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ABSTRACT 

An analytical technique for the calculation of linear decay 
coefficients in combustors with acoustic absorbers is presented. 
Tuned circumferential slot acoustic absorbers are designed for the 
first three transverse modes of oscillation and decay coefficients 
for these absorbers are found as a function of backing distance 
for seven different chamber configurations. The effectiveness of 
the absorbers for off design values of the combustion response 
and acoustic mode is also Investigated. ResultB indicate that for 
tuned absorbers the decay coefficient increases approximately as 
the cube of the backing distance. For most off design situations 
the absorber still provides a damping effect. However, if an 
absorber designed for some higher mode of oscillation is used to 
damp lower mode oscillations a driving effect is frequently found. 
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Nomenclature 

speed of sound 

absorber backing distance 

quantities defined after Equation (A-3) 

quantity defined after Equation (A-3) 

decay coefficient 

effective decay coefficient 

functions defined after Equation (A-A) 

Green's function satisfying Equation (A-2) 

modified Green's function defined after Equation (A-2) 

function defined after Equation (A-A) 

unit imaginary = /-T 

imaginary part of ( ) 

positive integer 

Bessel Function of the first kind of order m 
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positive integer 
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effective absorber aperture length 
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outward unit normal vector 
pressure 
velocity vector 
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chamber radius 

- radial length 

- position vector 
source position vector 

- real part of ( ) 

surface area of chamber 

- surface area of absorber slot 

surface area of nozzle exit plane 

- time 

- axial velocity 
chamber volume 

- radial velocity 

- circumferential velocity 

- absorber slot width 
axial coordinate 

absorption coefficient 
acoustic admittance 

- ratio of specific heats 

- acoustic eigenvalue for mode N 

- circumferential coordinate 

- imaginary part of complex frequency 

I 

- root of J (X. ) = 0 

m Jim 

- normalizing constant defined after Equation (A-2) 

- matrix discussed after Equation (A-4) 
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p - density 

T - sensitive time lag 

4> - perturbation velocity potential 

ip - normalizing constant defined after Equation (A-3) 

- acoustic eigenfunction for mode N 

U) - complex frequency of oscillation 

Subscripts 

A - absorber aperture quantity 

C - quantity evaluated on chamber cylindrical periphery 

i - injector quantity 

N - nozzle quantity 

o - mean value 

Superscripts 

- perturbation quantity or derivative with respect to 
argument 

* - dimensional quantity 

~ - chamber quantity with no absorber present 

- designates particular mode integers chosen 

(j) - approximation of order j 
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SUMMARY 

The effectiveness of a certain type of acoustic absorber as a damp- 
ing device in seven different rocket motors has been evaluated analyti- 
cally. 

The type of absorber considered in this study was an annular slot 
of variable backing distance and aperature length and width. The slot 
was assumed to be located immediately downstream of the injector. Tuned 
absorbers were designed for three different backing distances for both 
the first and third transverse modes of oscillation, and the decay rates 
for the chamber under consideration were then calculated. The effective- 
ness of each absorber design was also evaluated for three types of off 
design situations. The first of these was the case where the real part 
of the combustion response factor was larger than the design value. The 
second was the case where the real part of the combustion response factor 
was smaller than the design value. The third type of off design situa- 
tion was for a chamber mode of oscillation different from the design mode. 

Results of the study indicated a very strong dependence of the decay 
rate on the backing distance. It was also found that as either Mach 
Number or length to radius ratio for the chamber increased the combustors 
became relatively more unstable for a given backing distance. The calcu- 
lations indicated that a first transverse mode design absorber was more 
successful at damping third transverse mode oscillations than was a third 
transverse mode design at damping first transverse mode oscillations, 
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and that, in fact, third transverse mode designs could drive first trans- 
verse mode oscillations. The absorbers were, in general, less effective 
at off design combustion response values, though with sufficiently large 
backing distances the damping effect of the absorbers was often large 
even in these cases. 

Decay rates (log damping rates) were calculated using an integral 
equation iteration technique which calculates successively better approxi- 
mations for the velocity potential and complex frequency. The analytical 
technique does not restrict the magnitude of the mean chamber Mach Number 
and is designed to accommodate the discontinuous boundary conditions which 
arise when a partial length acoustic absorber of the type used in this 
analysis is present. Convergence of the numerical technique used was 
sufficiently rapid to recommend the future use of this analytical technique 
in the evaluation of acoustic absorber performance in future or existing 
rocket combustion chambers. 
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INTRODUCTION 

It is well known that acoustic absorbing devices can have a strong 
positive effect on the combustion stability of liquid rocket motors. 

It is therefore very desirable to have a method of evaluating analyti- 
cally the potential effect of such acoustic absorbers on the stability 
of rocket engines. 

Studies of the effects of full length acoustic liners on the 
stability of combustion chambers have been performed both for the simple 
case of a chamber with no mean flow, nozzle, or combustion zone (see 
Ref. 2) and for the more complex case of a chamber with finite mean flow, 
concentrated combustion, and short nozzle (see Ref. 1). More recently a 
method of evaluating the effects of partial length acoustic liners has 
been developed which also includes the effects of mean flow, a short 
nozzle and a concentrated combustion zone (see Ref. 4). This latter work 
is of the greatest practical interest since it is usual that acoustic 
absorbers line only a fraction of the chamber length in real engines. 

The results of this last mentioned study were given in terms of 
stability curves for the chamber using a combustion response factor 
which characterized the combustion process. For design and evaluation of 
rocket engine stability, however, it is often more useful to know the 
damping rate produced by a given absorber in a given engine for a given 
combustion response. It is necessary to modify the method mentioned 
above in order to give results in this form, and to include a frequency 
sensitive absorber impedance. It is the purpose of this report to pre- 
sent the modified theory and to show the results of applying it to seven 
rocket engine designs. 
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The type of absorber considered in this work was a simple circum- 
ferential slot of variable backing distance, aperture length, and 
aperture width. (See Fig. 1). The seven engines evaluated provided a 
rather large cross section of chamber design parameters such as mean 
Mach Number, length to radius ratio and ratio of specific heats. The 
influence of these parameters on the stability of chambers with acoustic 
absorbers of this type could therefore be determined. 

Decay rates are calculated using a modification of the method of 
analysis mentioned above and in Ref. 4. This modified analysis is an 
integral equation iteration technique in which the iterations are carried 
out numerically on a digital computer. The oscillations are assumed to 
be linear (small amplitude) , and to occur in a homentropic and irrota- 
tional source free field of calorically perfect gases. The combustion 
processes are represented using the sensitive time lag model. 

Results are given in terms of log damping rates for a given engine- 
absorber combination and combustion response (location on the neutral 
stability curve for the unlined chamber). Curves for only the first 
and third transverse primary modes of oscillations are presented here, 
although some calculations have also been performed for the second and 


fifth transverse modes. 
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THEORY 

Combustor Model 

Combustion chambers of cylindrical cross section terminated by a 
"short" nozzle and having all combustion concentrated at the injector 
face are the type considered in this report. The combustion zone mass 
generation response is assumed to be pressure dependent only and can be 
represented using the sensitive time lag model devised by Crocco (3) . 

An acoustic absorber is assumed to be located immediately downstream of 
the injector. This absorber is basically a radially oriented circum- 
ferential slot. Its geometry and location are shown in Figure 1. The 
rest of the cylindrical periphery of the chamber is taken to be non- 
absorbing. The gasdynamic field downstream of the combustion zone is 
assumed to be composed of a single component, calorically perfect com- 
bustion product gas. The flow is taken to be homentropic and irrotational. 

Using these assumptions the conservation equations for mass and 
momentum can be written down. (See for example Reference 1). These 
equations are then linearized using the small perturbation assumption. 

Thus p=l+p', p = 1 + p', a = 1 + a', u = M + u', v = v', w = w' 
where all state variables are non-dimens ionalized by their mean values, 
velocities are made non-dimensional by division with the mean sound speed 
and primes indicate perturbation quantities. The perturbed variables are 
represented as the product of a space dependent function and e^ 1 " where 
a) is the complex frequency (oo = Re(w) + iX) and t is time. Frequency is 
non-dimensionalized by multiplication with the ratio of the chamber radius 
to the mean sound speed, time is non-dimensionalized by division with the 
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same ratio, and space variables by division with the chamber radius. 

Because of the irrotationality assumption it is possible to reduce 
the perturbed conservation equations to a single partial differential 
equation in <j> , the space dependent part of the perturbation velocity 
potential. This equation and the expressions relating <p to the other 
perturbed dependent variables are given in Appendix A. 

Introduction of a Green's function for the chamber and the use of 
Green's Theorem make possible the transformation of this partial differ- 
ential equation and its associated boundary conditions into the two 
governing integrals equations which follow. 

♦ - n » *fff ( m2 0 + 2tH “H ) av o + /A ( ^(- + ) dS ° 

v s 

( 1 ) 



In the equations above the volume integrals are taken over the complete 
chamber volume and the surface integrals over all the boundary surfaces 
of the cavity (cylindrical walls, injector face and nozzle exit plane). 
The quantities 8, M and G^ which appear in the equations have the 

following meanings: 

8 specific acoustic admittance over which the integral 

V — I 

is being taken. For the nozzle exit plane 8 N = )• 

For the cylindrical solid walls 8 r = 0. For the 

absorber slot width, W , 8 = , where K is the 

A C YK 
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specific absorber impedance and is a function only 
of Re((o) when the absorber geometry is fixed. For 
the injector plane (combustion zone) 

6 i = M[n(i-e -1WT ) - . 

M mean chamber Mach Number. 

one of the acoustic eigenfunctions for the chamber 
with no mean flow or liner. 

ri N the eigenvalue associated with the eigenfunction 

G„(r/r ) the Green's function for the chamber. (See 
N o 

Appendix A for exact definition) 


Solution of the Equations 

An iterative technique is used to solve the integral equations for 
(jo and tp. First an initial functional form for <J> and an initial value 
for w are substituted into the integrals on the right hand side of the 
equations. The initial values used in this work were the velocity 
potential and frequency for the chamber with no liner. These quantities 
can be obtained using a separation of variables technique. Integration 
using these initial values then produces first approximations for to and 
<p . These quantities are in turn substituted into the integrals to yield 
second order approximations. This process continues until to and <}> are 
suitably invariant between successive approximations. Typically ten 
iterations yield results that vary less than one tenth of one percent. 

A more detailed description of the iteration process used is given in 
Appendix A. A description and print out of the computer program used 
in the iterative calculations is given in Appendix C. 
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RESULTS AND DISCUSSION 

Seven different engine designs were considered. Table 1 gives the 
necessary geometrical and mean operating parameters for these combustors. 
For each of the combustion chambers tuned absorbers (i.e. Im(K) = 0) were 
designed for two or more primary transverse modes of oscillation. All 
the absorbers were tuned for a combustion response value corresponding 
to the minimum point on the n,T neutral stability curve for the chamber 
without an absorber. An example of this point is shown in Figure (2) 
for Engine 1 as point A on the curve. 

The acoustic absorbers evaluated in this work were made to satisfy 
one of two design conditions. The majority of the absorbers were designed 
subject to the condition that the absorber aperture length be equal to 
the backing distance. Some of the absorbers, however, were designed in 
such a way as to maximize an "effective decay coefficient" for a given 
backing distance. This latter requirement in general vielded aperture 
lengths less than the backing distance. The "effective decay coefficient" 
was introduced in order to provide an approximate means of measuring the 
broad frequency band effectiveness of an absorber. Its exact definition 
and a discussion of its role in determining absorber designs in this 
work are presented in Appendix B. 

Once either one of these absorber design conditions was adopted for 
a particular engine design, absorber backing distance, and mode of 
oscillation, it was possible to calculate chamber frequencies and decay 
rates and absorber slot widths using the two governing integral equations, 
an expression for the impedance of the absorber (See Appendix B), and 
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the condition that the liner be tuned (Im(K) = 0). Tuned absorbers of 
this type were designed for several backing distances for a given chamber. 

If an absorber is designed so that it is tuned for a given chamber, 
a particular value of the combustion response, and a given primary mode 
of oscillation, then when either the combustion response value or the 
mode of oscillation is changed the absorber will not be tuned for the new 
conditions. For each absorber design studied at least two off design 
combustion response values and one off design mode were investigated in 
order to determine the general effectiveness of the absorber away from 
its tuning point. 

Tuned Absorber Results 

Curves of decay coefficient (imaginary part of the complex frequency) 
as a function of backing distance to chamber radius ratio for tuned first 
transverse mode absorbers for the seven chambers considered are shown in 
Figure (3) . The curves for each of the seven combustors are seen to 
approximate straight lines on the log-log plot. The slope of the lines 
is approximately three, indicating that the decay coefficient increases 
as the cube of the backing distance. All the absorbers used in Figure (3) 
were designed subject to the condition L a = B. The straight lines for 
the different engines, though approximately parallel, are in some cases 
displaced from each other a considerable distance. This means that the 
value of the decay coefficient for a given backing distance will be much 
lower in some engines than in others. For example, comparison of the 
lines representing Engines 3 and 4 shows that at a B/R value of 0.10 the 
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decay coefficient for Engine 3 is about 290 sec ^ while for Engine 4 it 
is only about 14 sec \ The difference in liner effectiveness is related 
to chamber design parameters such as Mach Number and length to radius 
ratio as will be shown shortly. 

The dependence of the decay coefficient on absorber slot width for 
the same engines, absorbers and mode as in Figure (3) is shown in Figure 

(4) . The slot widths shown for any engine are those that tune an absorber 
in that engine for the backing distances of Figure (3) . The absorber 
geometry in any engine for a particular decay coefficient is found using 
Figures (3) and (4) by simply picking the decay coefficient desired, 
noting the point of intersection of the line of constant decay coefficient 
with the desired engine line, and then reading off the corresponding 
values of B/R and W a /R. 

An attempt was made to correlate all the tuned absorber decay 
coefficient results for first transverse mode on a single plot. Figure 

(5) displays the results of this correlation when decay coefficient is 
plotted against the parameter 


a 



1 

Ml/ 3 



It is seen that all the absorbers evaluated fall fairly close to a line 
of slope approximately one. Since the plot is on a log-log scale, the 
scatter of the design points makes it clear that the attempted correlation 
is only roughly successful. It does, however, provide a means of obtain- 
ing an estimate of the backing distance required in order to obtain a 
given decay coefficient. It also indicates the effects of length to 
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radius ratio and Mach Number on the decay coefficient. Increasing 
either of these parameters decreases the decay coefficient for a fixed 
backing distance. The dependence is linear in the case of length to 
radius ratio and a 1/3 power dependence in the case of Mach Number. 
Because of the wide range of Mach Numbers and length to radius ratios 
examined, it is felt that some confidence can be placed in these depen- 
dencies. It should be remembered that the correlation of Figure (5) is 
valid strictly for tuned first transverse mode absorbers with L a = B. 

Third transverse mode tuned absorbers with three or more different 
backing distances were only designed for three of the engines studied, 
though at least one third transverse mode absorber was designed for all 
seven engines. Decay coefficients as functions of backing distance for 
the three engines evaluated most extensively are shown in Figure (6) . 

The three straight lines in the figure are those that result when the 
condition L a = B is used in the absorber design. The dashed lines in 
the figure represent absorbers designed subject to the maximum effective 
decay coefficient condition discussed in Appendix B. Until B/R was about 
0.09 it was impossible to satisfy this design requirement. (That is, the 
optimum effective decay coefficient occurred at L a values greater than B) 
Thereafter absorbers designed in such a way as to maximize the effective 
decay coefficient produced decay rates considerably smaller than those 
designed so that L a = B. An apparent trade off exists between these two 
types of absorbers at these larger backing distances in that an absorber 
designed for optimum broad band effectiveness produces a smaller decay 
coefficient under tuned conditions than does an absorber designed so that 
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L a = B. Comparison of Figure (6) with Figure (3) shows that, for the 
same backing distance, absorbers designed so that L a = B produce sub- 
stantially larger decay coefficients for third transverse mode oscilla- 
tions than they do for the first transverse mode. Figure (7) exhibits 
the dependence of the decay coefficient on the slot width for tuned third 
transverse modes. It can be used in conjunction with Figure (6) to 
determine absorber design dimensions, as was the case for first transverse 
mode absorbers. 

Off Design Results - Combustion Response 

The effectiveness of an absorber in damping oscillations occurring 
at combustion response (n,t) values different from design values was 
evaluated by determining decay coefficients for two points on the neutral 
stability curve for the given combustor with no absorber. One of these 
points was at a frequency higher than the design (minimum n) frequency, 

the other was at a frequency lower than the design frequency. Examples 

{ 

of these two points are shown in Figure (2) for the first transverse mode 
in Engine 1 as points B and C, respectively. Point A on the curve is 
the design point. Values of the real part of the frequency divided 

by the acoustic frequency for the first transverse mode, are shown as 
parametric points on the neutral stability curve. The effectiveness of 
the liner at the off design points varied between different chambers, 
different modes and different backing distances. For the example chosen, 
Engine 1, first transverse mode, backing distance B = 0.10, the decay 
coefficient at point A, the design point, was 123 sec 1 . At point B, 
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the high frequency point, it was -40 sec 1 . At point C, the low frequency 
point, it was 178 sec 1 . The neutral stability curves for this particular 
chamber with and without absorber are shown in Figure (8) . Two neutral 
stability curves are shown. The solid curve is for the chamber with no 
absorber; the dashed curve is the neutral stability limit for the chamber 
with the absorber. Areas inside neutral stability curves are linearly un- 
stable (growth occurs); areas outside are linearly stable (decay occurs). 
Since the dashed curve is predominantly inside the solid curve, the 
absorber provides an overall stabilizing effect for most of the frequency 
range of the first transverse mode. It is seen, however, that at fre- 
quencies above to / u) q = 1.04, the solid curve lies inside the dashed curve. 
For this region the addition of the liner is destabilizing. For example, 
point B, which is in this region, has a decay coefficient of -40 1 which 
indicates linear growth rather than decay. Without the absorber oscilla- 
tions at point B neither grow nor decay. The general trends in the re- 
sults of off design calculations of this type indicate that for frequencies 
above the design frequency the absorber is usually less effective and pro- 
duces smaller decay coefficients than at design. For very high frequencies 
close to the transition from a pure transverse mode to a combined trans- 
verse longitudinal mode it is possible for the absorber to be driving 
rather than damping, as in the example discussed above. For off design 
frequencies below the design frequency the absorber usually produced decay 
rates of the same order as those at design; in some cases larger than at 
design. Overall it was concluded that the stability improving effect of 
the absorber is quite dependent upon the value of the combustion response 
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and its location relative to the design value. In order to guarantee 
stability improvement for a given chamber, absorber, and mode combination, 
the entire range of combustion response factors of interest must be 
evaluated. 

Off Design Results - Higher Modes 

For all the first transverse mode absorbers designed, calculations 
were performed to determine the decay coefficients produced by them for 
the third transverse mode of oscillation. Some second transverse mode 
calculations of this type were also made. Typically the results of such 
calculations indicated that first transverse mode absorbers are quite 
effective in damping higher mode oscillations. Figure (9) shows the decay 
coefficients produced when a first transverse mode absorber was used to 
damp third transverse oscillations in three of the engines studied. Com- 
parison with Figure (3) shows that for small backing distances larger 
decay coefficients are produced for the third transverse mode than for 
the first transverse mode even though the absorber is tuned for first 
transverse mode oscillations. As backing distances become larger, how- 
ever, the fact that the slopes of the tuned first transverse mode lines 
(Figure (3)) are considerably larger than the slopes in Figure (9) 
ensures that their damping coefficients will eventually be larger. This 
effectiveness of the first transverse mode absorbers in damping third 
transverse mode absorbers is not too surprising in light of the fact that 
for the same backing distance tuned third transverse absorbers provide 
decay coefficients for the third transverse mode about an order of magni- 
tude larger than do tuned first transverse mode absorbers for the first 
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transverse mode (compare Figures (3) and (6)). 

Off Design Results - Lower Modes 

Calculations were also performed in order to determine how effective 
higher mode absorber designs were at damping lower mode oscillations. 
Results indicated not only that using higher mode designs for lower mode 
oscillations always produced rather small damping coefficients, but also 
that in many cases these coefficients actually were negative. That is, 
the absorber drove the lower mode oscillations . This type of situation is 
demonstrated in Figure (10). Figure (10a) shows decay coefficients as 
functions of B/R for the first transverse mode of oscillation with absorb- 
ers tuned for third transverse oscillations for three different engines. 

It can be seen that for backing distances less than about 0.15 the decay 
coefficients are negative. As the backing distance increases beyond 0.15 
all three engines eventually show some damping effect with the third trans- 
verse absorbers. The stability condition of Engine 2 with an absorber of 
backing distance 0.10 is shown in Figure (10b) using neutral stability 
limits. The dashed curve is a neutral stability limit for the first 
transverse mode using the absorber designed for the third transverse mode. 
The solid curve is the neutral stability limit for the engine with no 
absorber at all. It is seen that for the entire first transverse mode 
the solid curve is inside the dashed curve. This means that the third 
transverse mode absorber is destabilizing with respect to the first trans- 
verse mode for any value of combustion response. 

A higher mode absorber used to damp lower mode oscillations always 
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produced a negative imaginary part for the acoustic impedance of the 
absorber. Using a lower mode absorber to damp higher mode oscillations 
always produced a positive imaginary part for the acoustic impedance. 

The dependence of the non-dimensional decay coefficient on the imaginary 
part of the absorber impedance for a fixed value of the real part of the 
impedance is shown in Figure (10c) for the first transverse mode of oscil- 
lation in Engine 2. It is seen that maximum damping occurs at tuning, that 
is, where Im(K) = 0. Increasing the imaginary part of K decreases A until 
as A ■+ oo no damping effect is produced. Decreasing the imaginary part of 
K (Im(k) < 0) produces a driving effect that first increases and then de- 
creases to zero. It is to be noted that an extreme case has been presented; 
if Re (k) were smaller (more acoustic resistance) the entire curve in Figure 
(10c) would be shifted upwards and A would not necessarily ever be negative. 
The behavior of A, however, would always be asymmetrical with respect to 
Im(K) = 0 with values of Im(K) < 0 providing less damping than values of 
Im(K) >0. It is important to note that this type of asymmetrical behavior 

is not reflected using the absorption coefficient. As it is usually defined 

2 2 

(a = 4Re (K) / ([Re (K) + l] + [lm(K)J ) the absorption coefficient is abso- 
lutely symmetrical with respect to Im(K) = 0. The conclusion reached here 
is that the value of the absorption coefficient alone is not sufficient to 
determine the effectiveness of an absorber. For the same value of a the 
damping of a given absorber can be quite different depending on the sign 
of the imaginary part of the acoustic impedance. 
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SUMMARY OF CONCLUSIONS 

An analytical method for predicting decay coefficients in rocket 
combustors with partial length acoustic absorbers has been developed. 
Application of this method to seven different combustors for a particular 
type of absorber design led to the following conclusions. 

1. The technique was applicable and practicable for all combina- 
tions of absorbers, acoustic modes and engines tested. 

2. It was possible to correlate all the decay coefficients for 
the first transverse mode tuned absorbers using a correlation 
coefficient which involved the chamber length to radius ratio 
and mean Mach Number. From this correlation it was concluded 
that decay coefficients increased with the cube of the ab- 
sorber backing distance and decreased with Mach Number to the 
1/3 power and length to radius ratio to the first power. 

3. For values of the combustion response other than that value 
for which the absorber was tuned, the absorber was often less 
effective. In fact for some values of combustion response 
the absorber provided a driving rather than damping effect. 

4. When an absorber tuned for a particular mode was used to 
damp a higher acoustic mode it was quite effective. 

« 

5. When an absorber tuned for a given acoustic mode was used 
to damp a lower mode it was generally very ineffective and 
in many cases provided a driving rather than damping effect. 
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TABLE I 


Engine 

Number 

Mach 

Number 

Ratio of 
Specific Heats 

Length to 
Radius Ratio 

Radius of 
Chamber 

Speed of 
Sound 

1 

0.265865 

1.146 

1.984159 

9.07 in. 

5181. 7ft/sec. 

2 

0.304389 

1.148 

1.841964 

8.47 in . 

5248.1 

3 

0.392677 

1.149 

1.051433 

7.59 in. 

5281.6 

4 

0.154893 

1.343 

5.440007 

5.03 in. 

5569.6 

5 

0.081879 

1.3387 

3.732183 

4.44 in. 

5620.2 

6 

0.064638 

1.3835 

5.662040 

3.17 in. 

5512.3 

7 

0.191492 

1.348 

6.796209 

1.26 in. 

5740.6 
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APPENDIX A 
Basic Equations 

The perturbation equation describing the oscillatory flow in the 
combustion chamber written in terms of the perturbation velocity 
potential is 

V 2 cJ> + to 2 4> = M 2 f-f- + 2Mico 1^- (A-l) 

o z* d z 

This equation is derived from the conservation equations for the chamber 
in Reference (1). Perturbation pressure and velocity are related to <J> 
through the following relationships 

p' 2 — y (ico4> + M 

d z 

P-n 

Equation (A-l) can be transformed into the two integral equations 
Equations (1) and (2) through the introduction of the Green's Function 
satisfying the following equation 

V 2 G + w 2 G = 6 (r-r ) (A-2) 

o 

where $G • n = 0 on all boundaries and 6 is the delta function. The 
Green's function is represented as an expansion in tha orthonormal eigen- 
functions for the chamber with no combustion or mean flow and solid 
boundaries as 

o-EEZ n ^ (h 

ST m n w 2 - rin 

xmn 

where is one of the orthonormal eigenfunctions given by 
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Jlmn 


T /*\ \ mr z 

cos nc J (A- r) cos — — 
m £m L 


h 

£mn 


A 


and 


ilmn 


n 2 ir 2 

L 2 




The quantity A^ is a root of the equation 


L- 


0 and 


in X/ m 

^Jlmn 2 ' s determined by the condition J 'J “ £ ^£ mn ^ = 

The modified Green's function G^ which appears in Equations (1) 
and (2) is just the Green's function with one term in the series removed. 

A ^ ^ /N /V /\ 

That is», the term with il = H , m = m, n = n where H , m, n are arbitrary s 

is deleted from the series for G„. Thus 

N 


m n 


Sj. (r) S, (r o> 

nn X/titn 


s- ZEE- Cl - 6 (il ,t )6 (m,m)6 (n,n) 


.urn 


where <5(i,j) = 1 if i = j and 6(i,j) =0 if i ^ j . 

In Equation (1) the quantity can be interpreted as and in 

Equation (2) n„ is equivalent to nn /v/v • 

N X/inn 

In the solution of Equations (1) and (2) the velocity potential and 
frequency for the chamber with no liner but with mean flow and combustion 
are used as first approximations to (j) and u. These quantities are called 
respectively <p and u>. The form of $ can be found using the technique of 
separation of variables and is given by 


* = 


where 


B, = 


£ rr^- V x rm r)(eizBl +ce izB2 ) 

uM+ [w 2 M 2 + CM 2 -1) CX|- - w 2 )] 1 * 


(A-3) 


(1 - M 2 ) 
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+ 


y , 0 c G i| (u (j ” 1) )F2 W,0) (j-1) ) dS + $ J G N (w (j_1) )F 2 (<j> (j 1) -^)u) (j " 1) dS 


S-S 


/• 


+ /g c G N (co (j 1) )F 2 [ (4> (j-1) -|),W (j-1) ] dS 


+ M i:> )Fi[(4) (j 1) 


dV 


(A-4) 


co 




a) [v /ssyds + m fa* §f dv] 

L S V 


/en N F 2 [(cj) <j ' 1) 4),w (J L) ]d + My*Q N Fj[ (4) (;5 1} $),0) (j 1} ]dV 

S V 


where F 2 (4>,U3) = •Y(iu)<{> + M ||) 
Fi (4> ,00) = 21w || + M 0- 


(r Q ) 


H = ]C — (1-6 (m,m)6n,n) 

* T m a (a ) 2 (J ~ 1 ) -n§ f ) 


S-S is the surface of the two ends of the chamber. S is the surface 
c c 

of the cylindrical walls of the chamber. The expressions above are valid 
for j > 2. 

For j = 1 


<» - s + 


f e v 


<j> w = <t» + /8 G M (aj)F 2 (4,u) ds 
S 


03 


( “ ' “ 2 + / 


6 n N F 2 (o3) ds 
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B2 = <oM - [< 0 2 M + (M 2 -1) (X|- - w 2 )]* 5 
(1 - M 2 ) 

e iLBi [Bi + 6 N (YW + YMBi)] 

c = " ^IlbT [b 2 + B n (yw + y'mb 2 )] 

and 4> is determined by the condition 

/// » n r™ dv ■ 

V 

At the injector <j) must satisfy the contition ^<j> • n =■ 3^ p'» 

For a given combustion response this equation serves as an equation for 
the determination of to. This equation may be solved for to or, alternatively, 
Equations (1) and (2) may be solved for to after 8 C has been set equal to 
zero and <j> has been substituted for <$> . The latter method was used in this 
work . 

Substitution of $ and to into Equations (1) and (2) as the initial 

step in the iteration procedure described earlier eventually leads to the 

jth 

following two recussion formulas for the calculation of the approxima- 
tion for <J) and to 

<J> (J) + (w 2 [B j H N (w (;1 " 1) ,(o)F 2 ($,to (J ' 1) ) dS 

L S-S 

c 

+ M J* (w " 1) , to )Fi($,co (;J " 1 5 ) dV 
V 

+ i(to^ _1) - to) |y J G N (io) $ds + 2M G N (to) dv 

S-S V 
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After some quite extensive manipulation and itegration it is possible 
to rewrite the expression for in the following form 

4> (j) = $ + cos ae 2 £ C0S 

Z n 

y^ is a matrix of arbitrary dimensions. 

For the calculations in this report Z went from 1 to 3 and n from 
1 to 30. Larger matrices produced only slight increases in accuracy at 
great expense in computer time. . Using the 3 x 30 matrix fory^j^ resulted 
in computation times of between 1 and 2 seconds per iteration on a CDC 
6400 computer. 




Figure A-l Dependence of Decay Coefficient and Effective Decay Coefficient on Aperture Length. 
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APPENDIX B 
Absorber Design 

The following equation was used to calculate the impedance of the 
absorbers 

.. 0.1 


Y* 


+ i 


^Re ((jj) 


W 


£. 


'eff 2Re(co)B 2 


) ....(B-l) 


where £ ,, = L 
eff a 


r F“i 

+ 0.85 W I 1 - 0.7^2^- » and Re(w) is the real part of 


the chamber frequency* Simple Helmholtz resonator theory was used to 

obtain this equation. (See for example Reference (5)). For a tuned liner 

(0 . 1 
— - i — j . 

Clearly, for a fixed backing distance, B, and aperture length L , the 

o 

absorber can be tuned for a range of chamber frequencies by picking W Q 
ptoperly (so that lm(K) *» 0). Since the frequency and decay rate for the 
chamber are dependent not only on the absorber impedance but also on the 
aperture width, varying L for a given B will vary the chamber frequency 
and decay rate even though K remains fixed at its tuned value. For a 
given backing distance, then, frequencies and decay tatea (and absorber 
aperture width) for tuned absorbers can be calculated as a function of 
L . Note that in these calculations Re(w) appears in the imaginary part 
of R. This dependency must, of course, be included in the iterative 
solution method for Equations (1) and (2) mentioned earlier. 

For a given backing distance B there is some optimum aperture length 
L which will produce the maximum decay rate for the chamber. Unfortunately, 

Ei 

for the range of backing distances and chamber configurations examined in 


this report the optimum L usually had a value approximately an order of 

Gl 
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magnitude greater than the backing distance. Thus, for a backing distance 
of the order of one tenth the chamber radius optimum aperture lengths were 
of the order of the chamber radius. These aperture lengths were considered 
unrealistic. 

In order to develop a better design an "effective decay coefficient", 

D , was introduced defined by the following equation 
et i 

where D is the decay coefficient for the 
chamber, Ki is the impedance of the tuned absorber, K 2 is the impedance 
of the absorber at a frequency 1.5 times as large as the tuned frequency, 
and K 3 is the impedance at a frequency 2/3 as large as the tuned frequency. 
This effective decay coefficient serves as an approximate means of measur- 
ing the broad frequency band effectiveness of a given absorber design. 

For any absorber design the values of were always somewhat smaller 

than D. Moreover, optimum values of L a in the sense of maximizing 
occurred at substantially lower values of than did the optima in terms 
of maximum D. Still, the optimum L a for maximum was usually too 

large for practical purposes and the arbitrary design condition L = B was 
adopted in these cases. For third transverse and higher modes of oscilla- 
tion the optimum often occurred at a value less than B. When this 
happened, this smaller value of L was used in the absorber design. This 
latter type of behavior is shown in Figure (A-l), for the third trans- 
verse mode of oscillation in Engine 1 with a backing distance B «* 0.10. 
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APPENDIX C 
COMPUTER PROGRAMS 

This section describes how the computer programs are used to deter- 
mine the decay coefficient for a given set of input parameters which 
characterize the combustion chamber. The first program, COMBRES, cal- 
culates the neutral stability curves for the combustor. The second pro- 
gram, DECOEFF, uses the minimum interaction index which was determined in 
COMBRES to calculate the decay coefficient and absorber dimensions for a 
given backing distance. 

These programs were written in FORTRAN IV to run on a CDC 6400. 

Many portions of the programs and the nomenclature in the programs are 
the same. When the programs were written very little effort was made to 
optimize the programs with respect to running time. 

COMBRES 

This program calculates three neutral stability curves for a given 
combustor for the combined first radial * first, second, and third trans- 
verse acoustic modes, respectively, as a function of frequency. A com- 
bustion response factor is calculated for each frequency and is then con- 
verted into the corresponding interaction index and sensitive time lag 
values. The frequency is incremented until a minimum interaction index is 
determined. 

The input information required is: 

1. LENGNO - the problem identifier 

2. ALENGTH - the ratio of the length of the combustor to the radius 

of the combustor 


3. AMACH - mean flow Mach Number 
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4. GAMMA - ratio of the specific heats 

The output from this program includes the interaction index from 
which the minimum value can be determined. The combustion response factor 
corresponding to the minimum interaction index is then used as an input 
variable in DECOEFF. 

DECOEFF 

This program calculates the decay coefficient and absorber dimensions 
for a given backing distance of the absorber. This program uses the 
method of successive approximation to solve for the decay coefficient. A 
first guess at the frequency is calculated for the given combustion re- 
sponse factor. This frequency is used to calculate a new velocity potential 
and then a new frequency is calculated. This technique proceeds until the 
error between successive iterations is within the specified limits. The 
accuracy of the decay coefficient is limited only by the size of the co- 
efficient matrix used in the calculation of the velocity potential. 

The input parameters are: 

1. LENGNO - the problem identifier 

2. ALENGTH - the length to radius ratio for the combustor 

3. AMACH - mean flow Mach Number 

4. RADIUS - radius of the combustor in inches 

5. SOS - speed of sound in the combustor 

6. GAMMA - ratio of specific heats 

7. FREQ - the real part of the frequency ratio for the minimum 


interaction index. 



8. BETA! - combustion admittance corresponding to the minimum 
interaction index 

. 9. MHAT - transverse acoustic mode number 

10. NHAT - longitudinal acoustic mode number 

11. MATRIX - identifier used to determine whether or not the 

coefficient matrix is printed out 

12. B - ratio of the backing distance of the absorber to the 

radius of the combustor 

13. ERR0R1 - acceptable error in the real part of the frequency 

14. ERR0R2 - acceptable error in the imaginary part of the frequency 
15* LDEX - specifies the number Of terms in the radial direction 

to be calculated 

16. NDEX - specifies the number of terms in the longitudinal 
direction to be calculated 

The output from this program lists the input variables, the decay 
coefficient, and the absorber dimensions. 
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Sample Calculations 

This section is an explanation of how the two programs are used to 
design an acoustic absorber and calculate the resulting decay coefficient, 
subject to the condition that the aperature length equals the backing 
distance. The data cards for each sample calculation are shown at the end 
of the appropriate program listing. On the page immediately following 
each program is the corresponding listing of the output of the sample cal-* 
culation. For the sake of brevity only a partial listing of the output for 
COMBRES is given. 

The first step in designing an acoustic absorber and calculating a de- 
cay coefficient for the absorber designed is to calculate the combustion 
admittance as a function of frequency. This determines the neutral stability 
map for the combustor when no absorber is present and determines the inter- 
action index and sensitive time lag at the minimum of the stability curve. 
These neutral stability maps are calculated by COMBRES for the first, second, 
and third transverse acoustic modes. The form of the input for COMBRES is 
as followst 


Columns 

Variable 

Type 

1-10 

LENGNO 

Alphanumeric word 

11-20 

ALENGTH 

Real number 

21-30 

AMACH 

Real number 

31-40 

GAMMA 

Real number 


On the page following the one displaying the above input data is a 
partial listing of the output from this program. From the lower half of 
the page the minimum for the interaction index, N, is determined by 
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inspection and the corresponding frequency ratio is noted, W/WO.= 0.98. The 
frequency ratio is then used to locate the block of output containing the 
combustion admittance, upper half of page. In the output shown the frequency 
ratio is incremented by one-hundredth. This increment size may be changed 
by substituting the appropriate increment size on eleventh card from the 
bottom on page [F = F + CMPLX (increment , 0.0)]. It should be noted that 
there is an output block similar to that at the top of the page for each 
frequency ratio and that the output displayed is only for the first trans- 
verse acoustic mode. The output for the other transverse modes is similar 
to the above. 

Once the minimum interaction index is determined to sufficient accuracy* 
the second step is to run DECOEFF to obtain the decay coefficient with 
respect to the combustion admittance previously calculated. In the 
sample run of DECOEFF given here the minimum n used was for W/WO “ .982 
with a corresponding value of 6^ of 0.0761 + 0.063751. These values are 
obtained when COMBRES is used with the increments for W/WO being in 
thousandths rather than in hundredths as in the sample above. The 
aperature width of the acoustic absorber is calculated according to the 
assumption that the aperature length equals the backing distance. The 
form of the input for DECOEFF is shown immediately following the list of 
the program. The input for the program in on two data cards as follows: 

Card 1 


Columns 

Variable 

Type 

1-10 

LENGN0 

Alphanumeric word 

11-20 

ALENGTH 

Real number 

21-30 

AMACH 

Real number 
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31-40 

RADIUS 


Real number 

41-50 

SOS 


Real number 

51-60 

GAMMA 


Real number 

61-70 

W/WO 


Real number 

Card 2 



- 

Columns 

Variable 


Type 

1-21 

real part of BETAWI 

Real 

number with exponent 

22-42 . 

imaginary part of BETAWI 

Real 

number with exponent 

43-44 

MHAT 


Integer number 

45-46 

NHAT 


Integer number 

47-50 

MATRIX* 


Alphanumeric number 

51-55 

B 


Real number 

56-60 

ERR0R1 


Real number 

61-65 

ERROR 2 


Real number 

66-70 

LDEX 


Integer number 

71-75 

NDEX 


Integer number 

The last page of 

the computer listing shows 

the form 

of the output of 


DECOEFF. The output includes most of the input variables, the decay co- 
efficient, the effective decay coefficient, the absorber dimensions, and 
the location of the aperature opening with respect to the injector. 

*If the coefficient matrix is to be printed out, then MATRIX is 
inputed as YESB, where B is a blank. 
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NOMENCLATURE 



DECOEFF 

COMBRES 

to 

OMEGA 

OMEGA 

to 

W 

Not applicable 

M 

AMACH 

AMACH 

A 

GAMMA 

GAMMA 

n 

Function ETA 

Function ETA 


Function BESPRIM 

Function BESPRIM 

W 

Function BESSEL 

Function BESSEL 

A £mn 

Function FNORM 

Function FNORM 


Function PSI 

Function PSI 

B 1 

BETA1 

BETAI 

B 2 

BETA2 

BETA2 

C 

C 

C 

6 i 

BETAI 

BETAI 

6 c 

BETAC 

BETAC 

6 N 

BETAN 

BETAN 

y £n 

MU 

MU 

K 

AK 

AK 

L eff 

LEFF 

Not applicable 

l a 

AL or LA 

Not applicable 

W A 

WA 

Not applicable 


B 


B 


Not applicable 
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PROGRAM C OMBRES! OUTPUTS 01, TAPE6=<XJTPUT , 1 NPUT= 101, TAPFI5= 1 NPUT ) 
COMPLEX BETA1 , BETA2, BETAN, BETAC, BETA! , BETAWl, OMEGA, C, AK1 , 

1 AIC2, PS12, N { 1 50 ) , WW01150) 

COMPLEX G, AN, AK, A, A4, A5, BETA, DUM1 , DUM2, PS!, DUMMY, G1, 

| r 9 p- 

COMMON / BL<A / BETA 1 . BETA2, AMACH, GAMMA, OMEGA, ALENGTH, C, AK1 
1, AK2, BETAN, RLMDA, BETA!, BETAC, BETAWl, X, LHAT, MHAT, NHAT, 

1 PS 12 
< = 1 

READ (5,200) LENGNO , ALENGTH, AMACH, GAMMA 
25 GO TO (20,21 ,22,23) ,K 

20 MHAT = 1 : P = CMPLX (0.90, 0.0 ) : GO TO 

21 MHAT = 2 : F = CMPLX (1.50, 0.0 ) : GO TO 24 

22 MHAT = 3 : F = CMPLX (2. 00, 0.0 ) : GO TO 24 

23 CALL EXIT 

24 K = < + 1 
1 = 1 

X = 0.0 

BETAC = CMPLX (0. 0,0.0) 

7 CONTINUE 
WWO(l) = F 

A = CMPLX ( 1.0 ,0.0 ) 

G = CMPLX ( 0.9166, 0.0 ) 

BETAN = AMACH* (G-1. /GAMMA) 

LHAT = 1 
NHAT=0 


RLMDA = BESPRIM ( MHAT + 1 , LHAT ) 

NHAI = NHAT + 1 
OMEGA=F* 1.84 11 8378 

BETA=CSQRT IOMEGA ’AMACH’ OMEGA ‘AMACH+ (AMACH' AMACH- 1.)* (RLMDA* *2 - 
IOMEGA ’OMEGA) ) / ( 1 . -AMACH ’AMACH) 

BETA1 = (OMEGA' AMACH) / ( 1 . -AMACH’ AMACH) +BETA 

BETA2= (BET A 1 ; -2. ’BETA 

A4 = CMPLX ( -AIMAG(BETAI) , REAL (BETA1 ) ) 

D'JMI = CMPLX! -A! MAG (BETA2) , REAL (BETA2) ) 

C = - CEXP (ALETJGTH’A4) / CEXP(ALENGTH’DUMI) 

C = C ’ (BETA1 + BETAN ’ GAMMA ’ (OMEGA + AMACH ’ BETA 1 )) / { 

1 BETA2 ♦ BETAN * GAMMA ' (OMEGA + AMACH ’ BETA2 J) 

PS i 2 = PS! ! DUWY ) "2 

AK1=!OMEGA+ AMACH ’BETA1+C* (QMEGA+ AMACH*BETA2) ) 'GAMMA 
AK1= CMPLX i -A 1 MAG (AK1), REAL! AK1 ) ) 

D'JMI = CMPLX (-AIMAG(BETAI), REAL (BETA1)) 

D’JMI = D'JMI ’ALENGTH 
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CJM2 = CMPLX (-AIMAG(BETA2), REAL (BETA2J) 

DUM2 = DUM2* ALENGTH 

AK2 = ( (OMEGA+AMACH'BETAI ) *CEXP(DUM1 )+ (0MEGA+AMACH*BETA2) *C*CEXP( 
1DUM2) ) 'GAMMA 

AK2 = CMPLX (-A1MAG (AK2) f REAL ( AK2) ) 

CALL CALCA4( A4, A5 ) 

BET AW 1 = (OMEGA 1 *2 - ETA (LHAT, MHAT, NHA1 .ALENGTH, RLMDA)- (AMACH*G2( 
1 NHAT ) /EPS 1 ZN (LHAT , MHAT , NHAT , ALENGTH, RLMDA ) +BETAN* AK2/EPS 1 ZTJ (LHAT , 
1 MHAT, NHAT, ALENGTH, RLTCA) * (-1 . )** (NHAT+2) ) /PS I (DUMMY) )/A4 
AN = 1. / GAMMA - BETAWl / AMACH 
IF ( REAL ( AN ) .LT. 0.00 ) GO TO 1 
N ( 1 ) = AN 
WRITE (6, 100)LENGNO 
WRITE(6,116) AMACH 
WRITE(6, 1 17) GAMMA 
WRITE (6, 118) OICGA 
WRITE (6, 119) ALENGTH 
WRITE (6, 101) LHAT, MHAT, MHAT 
WRITE(6, 1 15) X 
WRITE (6, 113) F 
WRITE (6,. 114) A 
WRITE (6, 102) BETAN 
WRITE (6, 120) G 
WRITE (6, 121) BETAC 
WRITE (6, 108) 

WRITE (6, 105) BETAWl 
WRITE (6, 103) AN 

BETAWl = ( BETA1 + C ' BETA2 ) / ( GAMMA * (OMEGA * ( 1 . + C ) + 

• 1 AMACH * ( BETA1 + C ' BETA2 ) ) ) 

WRITE (6. 1 12) BETAWl 

AN = 1 . / GAMMA - BETAWl / AMACH 

WRITE(6, 103) AN 

1 = 1 + 1 : F = F, + CMPLX ( 0.01 , 0.0 ) 

IF tREAL (AN) .LT.1.8 .AND. A1MAG (AN) .LT. 1.0 ) GO TO 7 
500 CALL NTAU(AMACH, REAL (AK),N,WWO, M) 

CALL EXIT 
. GO TO 25 

1 F = F + CMPLX (0. 01 , 0.0 ) 

GO TO 7 

100 FORMAT Cl THIS IS ENGINE NUMBER ’A10) 

101 FORMAT CO THE PRIMARY MODE ASSUMED IS LHAT = \!2,* MHAT = *,(2,‘ 
INHAT = *,!2) 

102 FORMAT CO BETAN = ’,2021.14) 
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103 FORMAT (* + ’ , 70X, * N = \2G21 . 14, ' 1 ') 

104 FORMAT!* \ 13,* * , 1 0G13.6) 

105 FORMAT (*0 BET AW 1 = \2G21.14) 

106 FORMAT (’0 N J EQUAL \ 12, * L EQUAL Ml,' TO M2) 

108 FORMAT ('+' , 70X, ' K = INFINITY’) 

109 FORMAT (*0 BETA) (*,12, ’) =‘,2G21.14) 

112 FORMAT (’0 BET AW ID = *,2G21 . 14) 

113 FORMAT l* 0 W/W0 = \2G21.14) 

114 FORMAT ( '+* , 70X, ' A = ’.2G21.14) 

115 FORMAT l’+\70X,* THE LENGTH OF THE LINER IS ’.G21.14) 

116 FORMAT ( 'O THE MACH NUMBER IS * , G2T . 14) 

117 FORMAT (*+’,70X,* THE RATIO OF SPECIFIC HEATS IS \G21.14) 

118 FORMAT ( *0 THE FREQUENCY IS ', 2G21.14) 

119 FORMAT l' + ',70X,' THE LENGTH OF THE COMBUSTOR IS \G21.14) 

120 FORMAT (’ + ’,70X,* G = ' , 2G2U4) 

121 FORMAT (*0 BETAC = *, 2G21.14) . 

200 FORMAT ( A 1 0 , 3F 10.0) 

END 

FUNCTION BESPR1M (M, L) 

DIMENSION A (10,5) 

C"" THESE are the ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF OR 
C”” SET EQUAL TO ZERO 

A(1 ,1) = 0.00000000 
A (2 .1) = 3.83170597 
A (3 ,1) = 7.01558667 
A (4 ,1) = 10.17346814 
A (5 ,1) = 13.32369194 
A (6 ,1) = 16.47063005 
A (7 ,1) = 19.61585851 
A (8 ,1) = 22.76008438 
A (9 ,1) = 25.90567209 
A ( 1 0, 1 ) = 29.04682853 

C*“* THESE are the ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF OR 
C"" SET EQUAL TO ZERO 

Ad ,21 = 1.84118378 
. AI2 ,21 = 5.33144277 
M3 ,21 = 8.53631637 
AI4 ,2) = 11.70600490 
AI5 ,2) = 14.86358863 
A (6 ,21 = 18.01552786 
A 17 , 21 = 21.16436986 
AI8 ,21 = 24.31132686 
A (9 ,21 = 27.45705057 
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AM 0,2) = 30.60192297 

C“” THESE are THE ROOTS or THE DERIVATIVE OF THE BESSEL FUNCTION OF OR 
C”*’ SET EQUAL TO ZERO 

All ,3) = 3.05423693 
A 12 ,3) = 6.70613319 
. A (3 , 3) = 9.96946782 
A (4 ,3) = 13.17037086 
A 15 ,3) = 16.34752232 
A (6 ,3) = 19.51291278 
A (7 ,3) = 22.67158177 
A (8 ,3) = 25.82603714 
A (9 ,3) = 28.97767277 
A (10,3) = 32.12732702 

C * * ‘ ' THESE ARE THE ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF OR 
C*“* SET EQUAL TO ZERO 

A(1 ,4) = 4.20118894 
A (2 , 4) = 8.01523660 
A (3 ,4) = 11.34592431 
A (4 ,4) = 14.58584829 
A (5 ,4) = 17.78874787 
A (6 ,4) = 20.97247694 
A (7 ,4) = 24.14489743 
A (8 ,4) = 27.31005793 
A (9 ,4) = 30.47026881 
A (1 0,4) = 33.62694918 

C*“* THESE ARE the ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF OR 
C**” SET EQUAL TO ZERO 

Ad ,5) = 5.31755313 
A (2 ,5) = 9.28239629 
A (3 ,5) = 12.68190844 
A 14 ,5) = 15.96410704 
A (5 ,5) = 19.19602880 
A (6 ,5) = 22.40103227 
A (7 ,5) = 25.58975968 
A (8 ,5) = 28.76783622 
. A (9- ,5) .= 31.93853934 
Ad 0,5) = 35.10391668 
1 BESPR'iM = A (l, M> 

RETURN 
ENTRY BESSEL 

C““ THESE ARE^THE BESSEL NUMBERS OF ORDER ZERO FOR JHE ZEROS OF THE B 
C*’** FUNCTION 

AM ,1) = 1.9Q00000C 
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A (2 ,1) = -0.4027588095 
A (3 ,1) = 0.301128303 
A (4 , n = -0.249704877 
A (5 ,1) = 0.218359407 
A (6 ,1) -0.19645371 

A (7 , V) = 0.180063375 
A (8 ,1) = -0.167184600 
A (9 ,1) = 0.156724985 
A ( 1 0, 1 ) = -0.148011108 

C*’“ THESE ARE THE BESSEL NUMBERS QE ORDER ONE TOR THE. ZEROS OF THE B 

C’“* FUNCTION 

A ( 1 ,2) = 0.5818649368 
A (2 ,2) = -0.3461258542 
A (3 ,2) = 0.2732981131 
A (4 ,2) = -0.233304416 
A (5 ,2) = 0.207012651 
A 16 ,2) = -0.188017488 
A (7 ,2) = 0.173459050 
A(8 ,2) = -0.161838211 
A (9 ,2) = 0.152282069 
A (10,2) = -0.144242905 

c ‘ * ' THESE are THE BESSEL NUMBERS OF ORDER TWO FORvTHE ZEROS OF THE B 

C**'* FUNCTION 

A(1 ,3) = 0.4864961885 
A (2 ,3) = -0.3135283099 
A (3 ,3) = 0.2547441235 
A (4 ,3) = -0.220881581 
A (5 ,3) = 0.197937434 
A (6 ,3) = -0.181010000 
A (7 ,3) = 0.167835534 
A (8 ,3) = -0.157195167 
A (9 ,3) =0.1 48363778 
A( 10,3) = -0.140878333 

c"" these are the bessel mjmbers of order three for the zeros of the b 

C"V -FUNCTION 

A(1 ,4) = 0.4343942763 
A (2 ,4) = -0.2911584413 
A (3 ,4) = 0.240738175 
A (4 ,4) = -0.210965204 
A (5 ,4) = 0.190419022 
A (6 ,4) = -0.175048405 
A (7 ,4) = 0.162954965 
A (8 ,4) = -0.153102409 
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A (9 ,45 =0.1 44866574 
AC 10,4) « -0.137844513 

C ’ * * * THESE ARE THE 8ESSEL NUMBERS OF ORDER FOUR FOR THE ZEROS OF THE B 
c”*’ function 

A(1 ,5) = 0.3906514545 
A (2 ,5) = -0.2743809949 
A (3 ,5) = 0.229590468 
A (4 ,5) = -0.202763849 
A (5 ,5) = 0.184029896 
A (6 ,5) = -0.169878516 
A (7 ,5) = 0.158655372 
A 18 ,5) = -0.149451156 
A (9 ,5) = 0.141714307 
A (10,5) = -0.135086328 
GO TO 1 
END 

SUBROUTINE CALCA4( A4, A5 ) 

COMPLEX BETA1 , BETA2, BETAN, BETAC, BET A 1 , BETAWl , OMEGA, C , AK 1 , 

1 AK2 , PS 1 2 

COMPLEX A4, A5. PS I . DUMMY, Cl, G2 

COMMON / BLKA / BETA1 , BETA2, AMACH, GAMMA, OMEGA, ALENGTH, C, AK1 
1, A<2, BETAN, RLMDA, BETA], BETAC, BETAWl, X, LHAT, MHAT, NHAT, 

1 PS 1 2 

A4 = A< 1 /EPS 1ZNI LHAT, MHAT, NHAT, ALENGTH, RLMDA) /PS 1 (DUMMY) 

A5=BETAC*G1 (NHAT) /EPS I KAT (LHAT, MHAT, NHAT, ALENGTH, RLM)A)/PS1 (DUMMY) 

return 

EM) 

FUNCTION EPS I KAT (L.M.N, ALENGTH, RLMDA) 

DUM1 = ALENGTH 

IF ( N .NE. Oi-'DUMI = DUMI/2. 

IF ( M .EQ. 0) GOTO 1 

EPSIKAT = (0.5-M *M /(2. *RlM)A*RLMDA) ) *DUM1 
RETURN 

1 EPSIKAT = DUMI/2. 

RETURN 

ENTRY EPS1ZN 

IF ( N .EQ. 0 ) GO TO 2 

EPSIKAT= ALENGTH / 2. 

RETURN 

2 EPS1KAT= ALENGTH 
RETURN 

ENTRY ETA 

EPSIKAT = (N-DMN-D ’ 3.1415926 * 3; 1415926 / (ALENGTH * ALENGTH) 
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1+RLHDA V RLMDA 
RETURN 
END 

COMPLEX FUNCTION G1 ( N ) '' 

COMPLEX BETA1, BETA2, BETAN, BETAC, BETA I , BET AW I , OMEGA, C, AK1 , 

1 A K2, PS! 2 

COMPLEX DUM1, GDUM, BETA, GBUM, DUM 

COMMON / BLKA / BETA1, BETA2, AMACH, GAMMA, OMEGA, ALENGTH, C, AK1 
1, AK2, BETAN, RLMDA, BETA], BETAC, BETAWl, X, LHAT, MHAT, NHAT, 

1 PSI2 

GDUMfA, DUM! .BETA, OMEGA, AMACH, X) = (OMEGA+AMACH'BETA) / (A"2~BETA"2 
1 ) * (CEXP(DUM1 * X ) * !DUM1 'COS(A'X) +A*S1N(A'X) ) -DUM1 5 
GBUM (BETA , DUM! , N, AMACH, OMEGA , ALENGTH) = - (AMACH'BETA'BETA+2. ’OMEGA 
1 'BETA) / ( ( N ' 3.1415926 / ALENGTH ) ”2 - BETA ”2 ) ’(CEXP(DUM)' 

1 ALENGTH) ' DUM! * (-1.)"(N+2) - DIM! ) 

DUM (DUM1, ALENGTH, NHAT) = CEXP (DUM1* ALENGTH) * (-1 .) ” (NHAT+2M . 

IF ( X .EQ. 0.0 ) GO TO 1 
A = N'3. 1415926/ALENGTH 
BETA = BETA! 

DUM1 = CMPLXI-A! MAG (BETA!), REAL (BETA!)) 

G1 = GDUM(A,DUM1, BETA, OMEGA, AMACH, X) 

IF ( N .EQ. 0 .AM). CABS ( BETA2) .LE. 1.E-10 ) GO TO 2 
BETA = BETA2 

DUM1 = CMPLX (-A1MAG (BET A2) .REAL (BETA2)) 

G! = GAMMA ' (G1 + C * GDUM ( A, DIM! , BETA,, OMEGA, AMACH, X) ) 

G1 = CMPLX ( -AIMAG(GI), REAL (GD) 

RETURN 

1 G1 = CMPLX(0. 0,0.0) 

RETURN 

2 BETA = C ' OMEGA ' X 

G! = GAMMA ' ( G1 + BETA ) 

G1 = CMPLX ( -AIMAG(GI), REAL (G 1 ) ) 

RETURN j 

ENTRY G2 
BETA = BETA1 

DUM! = CMPLX (-AIMAG(BETAI), REAL (BETA!)) 

G1 = GBUMiBETA, DUM! ,N, AMACH, OMEGA, ALENGTH) 

IF ( N .EQ. 0 .AND. CABS ( BETA2) .LE. 1.E-10 ) RETURN 
BETA = BETA2 

D'JM! = CMPLX (“AIMAG(BETA2), REAL (BETA2)) 

. Cl = G! + C * GBUM (BETA, DUM1.N, AMACH, OMEGA .ALENGTH) 

return 

ENTRY PS! 
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DUM1 = CMPLXl -AIMAG (BETA 1 ) , REAL (BET A 1 ) ) 

G1 = DUM1 *DUM(DUM1 , ALENGTH, NHAT) /(NHAT’ *2*3. 1415926* *2/ ALENGTH* *2 
1 -BETA1 * *2) 

IF t NHAT .EQ. 0 .AND. CABS ( BETA2) .LE. 1.E-10 ) GO TO 3 
DUM1 = CMPLXl -AIMAG (BETA 2),- REAL (BETA2) ) 

G1 = ( G1 +C*DUM1 ‘DUMtDUMI .ALENGTH, NHAT) / (NHAT* *2*3. 1415926* *2/ 

1 ALENGTH’ *2- BETA2 **2)1 /EPS I ZN (LHAT , MHAT , NHAT , ALENGTH, RLMD A) 

RETURN 

3 G1 = ( G 1 + C ’ ALENGTH ) / EPS 1 ZN (LHAT , MHAT , NHAT , ALENGTH, RLMD A , 

RETURN 
EW 

SUBROUTINE NTAU(MACH,K, AN, WWO, I) 

COMPLEX AN ( 1 ) , WWOd) 

REAL N, MACH, K, LINER 
WRITE (6, 102) 

P1E2' = 2. *3.14159265358979323846 
DO 7 J = 1,1 
ANR = REAL (AN (J)) 

AN1 = AIMAG (AN(J)) 

0 = REAL (WWO (J) ) 

OMEGA = 0 * 1.84118378 
N = ( ANR* *2 + AN! * *2 ) / 2. / ANR 
IF ( ANI .LT. 0.0 ) 5, 6 

5 TAU = (P I E2-AC0S ( 1 . -ANR/N) ) /'OMEGA . 

GO TO 7 

6 TAU = ACOS ( 1 . -ANR / N ) / OMEGA 

7 WRITE (6. 101) MACH, 0, OMEGA, N, TAU 

101 FORMAT (*0\5X,F5.3,5X,F5.2 ,5X,G18. 1 1 ,2(5X,G21 . 14) ) 

102 FORMAT (* 1 *,6X, *MACH*,6X, *W/W0*,9X, ‘OMEGA*, 20X, *N*,24X, ’TAU*) 

RETURN 

END 


THE FOLLOWING IS AN EXAMPLE OF THE PUNCHED CARD INPUT. 
ONE 1.984159 0.265865 1.146 


\ 
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THIS IS ENGINE NJMBER ONE 

THE MACH MJMBER IS .26586500000000 

THE RATIO OF SPECIFIC HEATS IS 1 . 1460000000000 

THE FREQUENCY IS 1.8043601044000 0. 

THE LENGTH OF THE COMBUSTOR IS 1.9841590000000 

THE PRIMARY MODE ASSUMED IS LHAT = 1 MHAT = 1 NHAT = 0 

THE LENGTH OF THE LINER IS 0. 

W/W0 = .97999999999998 0. 

A = 1.0000000000000 0. 

BETAN = 1.16979672024435E-02 0. 

G = .91660000000000 0. 

BETAC = 0. 0. 

K = INFINITY 

BETAWI = 8. 2532350 141 4754E-02 7. 0091 2836202505E-02 
N = .56217080719944 -.26363486589153 1 

BETAW1D = 8. 25323501 41 471 9E-02 7. 00912836202479E-02 
N = .56217080719946 -.26363486589152 I 


MACH 

w/wo 

OMEGA 

N 

TAU 

.266 

.90 

1.66 

20.526968 

3.7584663 

.266 

.91 

1.68 

7.4276621 

3.6665508 

.266 

.92 

1.69 

3.7860679 

3.5616498 

.266 

.93 

1.71 

2.1566187 

3.4374212 

.266 

.94 

1.73 

1.2912086 

3.2845350 

.266 

.95 

1.75 

.80152975 

3.0904435 

.266 

.96 

1.77 

.52708698 

2.8429163 

.266 

.97 

1.79 

.38891625 

2.5435533 

.266 

.98 1 

1.80 

.34290233 

2.2271651 

.266 

.99 

1.82 

.36160715 

1.9513585 

.266 

1.00 

1.84 

.42616618 

1.7525075 

.266 

1.01 

1.86 

.52232427 

1.6324354 

.266 

1.02 

1.88 

.63840129 

1.5780414 

.266 

1.03 

1.90 

. 76426740 

1.5756145 

.266 

1.04 

1.91 

.89090300 

1.6137839 

.266 

1.05 

1.93 

1.0103177 

1.6820456 
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PROGRAM DECQEFF ( I NPUT =101 , OUTPUT = 1 0 1 , T APE5= I NPUT, , T APE6=0UTPUT ) 
COMPLEX BETA1, BETA2, BETAN, BETAC, BETA!' BETAWl , OMEGA, C, AK1 , 

1 AK2, PS12, MU, W, F, 1 , SOME , SOME 1 , SO 1 , SSG 1 , Y , Z 
COMPLEX G, AN, A<, A, A4, A5, BETA, DUM1 , DUM2, PSi, DUMMY, G1, 

1 G2, Ml ,UHAT 

COMMON / BLKA / BETA 1 , BETA2, AMACH, GAMMA, OMEGA, ALENGTH, C, AK1 
1, AK2, BETAN, RLMDA, BETA 1 , BETAC, BETAWl, I, LHAT, MHAT, NHAT, 

1 PS12 , PI, X2, XI, SOME, SOME1 
COMMON / BLKB / MU! 3,30,29) , W(30) 

COMMON / BLK1 / E,RHOCHAM,RHOVOL,PRESURE,RADIUS,SOS. 

COMMON Cll (30,30) , CIO (30,30) , ERROR! , ERR0R2 , LDEX , NDEX 
8 CONTINUE 

1 = CMPLXIO.O ,1.0 ) 

PI = 5.14159265358974 
A = CMPLX ( 1.0 , 0.0 ) 

G = CMPLX ( 0.9166, 0.0 ) 

READ (5,200) LENGNO , ALENG TH , AMACH , R AD 1 US , SOS , G AMMA , FREQ , BET A I , MHA T , 
INHAT , MATR 1 X , B, ERROR 1 , ERROR2 , LDEX , NDEX 
IF ( EOF ,5 ) 9 , 10 
10 RAD * RADIUS 

AN = 1 . / GAMMA - BETA 1 / AMACH 
RADIUS = RADIUS/12.0 
F = CMPLX ( FREQ , 0.0 ) 

4 W = REAL (F)* 1.841 18378 
AP = 0.595*0. 595/2. /B 
BP = 0.85-0.5/ (W*B) ”2 

CALL CUB1C(-BP*BP/AP,-2.'B*BP/AP,-B*B/AP,B,WA) 

CALL ABSORB ( AN,B, WA , AL , AK , UHAT , GAMMA) 

FREQ1 = REAL(F) . 

XI = ( B-WA/2.) 

X2 = XI + WA i 

BETAC = 1,/AK /GAMMA 

BETAN = AMACH * ( G - 1 . / GAMMA ) 

LHAT = 1 

RLMDA = BESPR1M ( MHAT + 1 , LHAT ) 

7 CONTINUE 
NHA1 = NHAT + 1 
F= C^LX ( FREQ , 0.0 ) 

OMEGA = CMPLX ( 1.84118378 , 0.0 ) * F 

DO 2 1K= 1 , 500 

8ETA=CSQRT (OMEGA* AMACH *OMEG A 'AMACI-H- (AMACH* AMACH- 1 . ) * (RLfCA **2 - 
10f€GA*0f€GA) )/( 1. -AMACH* AMACH) 

BETA1 * (ONtGA* AMACH) / ( 1 . -AMACH* AMACH) +BETA 
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BETA2=(BETA1)-2.*BETA 
C = - CEXP (ALENGTH* 1 ' (BETA1- BETA2U 1 

C = C * (BETA1 + BETAN * GAMMA * (OMEGA + AMACH * BETA! )) / ( 

1 BETA2 + BETAN * GAMMA * (OMEGA + AMACH * BETA2 ) ) 

PS 12 = PS1 ( DUMMY ) * '2 

AK1= (OMEGA+ AMACH ‘BETA1+C* (OMEGA+ AMACH ' BETA2) ) 'GAMMA* I 

DUM1 = BETA! * I ‘ALENGTH 
DUM2 = BETA2* I ‘ALENGTH 

AK2 = ( (OMEGA+AMACH*BETA1 ) ‘CEXP (DUM1 ) + (0MEGA+AMACH'BETA2) *C‘CEXP( 
1DUM2)) ‘GAMMA*! 

W 1 = ( ( BET AN ’ AK2 * C OS ( NHA T '' 3 . 1 4 1 59265358974 } +BET A 1 ’ AK1 + AMACH * G2 ( NHA 
IT)) /PS 1 (NHAT) /EPS I ZN(L,M,MHAT, ALENGTH, RLMDA)+ETA(L,M,MHA1, ALENGTH, 
1RLMDA) ) 

Wl = CSQRT(Wl) 

IF ( REAL ( W1 ) .LT. 0.0 ) Wl = - Wi 
PS 12 = PS1 ( DUMMY ) 

IF ( CABS (Wl -OMEGA). LT.1.E-07) 3 , 2 

2 OMEGA = Wl 
WR1TEI6.125) 

CALL EXIT 

3 OMEGA = Wl 

F = Wl / 1.84118378 
CALL CALCA4I A4, A5 ) 

BETAW1 = (0MEGA ; ‘2 - ETA(LHAT,MHAT,NHA|, ALENGTH, RLMDA)-(AMACH‘G2( 

1 NHAT ) /EPS I ZN (LHAT , MHAT , NHAT , ALENGTH, RLMDA ) +BETAN* AK2/EPS I ZN ( LHAT , 
1MHAT.NHAT, ALENGTH, RLMD A) * (-1 . ) “ (NHAT+2) )/PS! (DUMMY) )/A4 
WRITE (6, 100) LENGNO 
WRITE (8, 1 16) AMACH 
WRI TE (6, 1 1 7) GAMMA 
WRITE (6, 118) OMEGA 
WRITE (6, 1 19) ALENGTH 
WR1TE(6,10D LHAT, MHAT, WAT 
WRITE (6, 1 15) XI , X2 
WRITE (6, 1 13) F 
WRITE (6, 1 14) A 
WRITE (6, 102) BETAN 
WRITE (6, 120) G 
WRITE (6, 121 ) BETAC 
WRITE (6, 108) AK 
WRITE (6, 105) BET AW 1 
AN = 1 . / GAMMA - BETAWl / AMACH 
WRITE (6, 103) AN 

BETAWl = ( BETA1 + C ' BETA2 ) / ( GAMMA ‘ (OMEGA * ( 1 . + C ) + 
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1 AMACH * ( BETA! + C * BETA2 ) ) ) 

WR 1 TE (6, 1 12) BETAWl 

AN = 1. / GAMMA - BETAWl / AMACH 

WRITE (6, 103) AN 

AN = 1. / GAMMA - BETA! / AMACH 

WRITE (6, 122) BETA 1 

WRITE (6, 103) AN 

WRITE (6, 126) OMEGA , F 

WRITE (6, 130) B , WA , UHAT , AL 

WRITE (6, 131) FREQ1 

IF ( ABSIAIMAG (BETAC) ) .GT. 1.E-05 ) GO TO 1 

WI= (BETAC* ( G 1 (X2,NHAT)-G1 (XI ,NHAT) ) /PS 1 2 /EPS1KAT (LHAT.MHAT.NHAT, 
1ALENGTH.RLM0A) )+0MEGA*’2 
Z=W1 
NN=0 

WI = CSQRT(WI) 

IF (REAL (WI ) .LE.0.0JWI = -WI 
Y=WI 

SOME=SG1 (X2.NHAT) - SGI (XI ,NHAT) 

SOME 1 =SSG 1 (ALENGTH.NHAT) 

SOME= ( ( (BETA I * ( 1 . +C) + BETAN* M . ) ”NHAT’ (CEXPIDUM1 ) + C’CEXPl 
1DUM2) ) ) * 1 ’GAMMA +2. ’ 1 ’AMACH’SOMEI ) /EPS I ZN (LHAT.MHAT.NHAT, ALENGTH, 
2RLMDA) + I* GAMMA * BET AC ’ SOME/EPS I KAT (LHAT , MHAT , NHAT , ALENGTH , RLMDA ) 
3) /PS 12 
14 CONTINUE 

WI=Z + (Y-QMEGA) ‘SOME 
WKSQRT(WI) 

IF (REAL (WI ) .LE.0.0) WI=-W1 
IF (CABS (WI - Y).LT.1.E-07.0R.NN.GT.10) 12,13 
13 Y=WI 
NN=NN +1 
GOTO 14 
12 W(1)=W1 

F = WI / 1.84118378 
WR 1 TE ( 6 , 109) Wl 
WRITE (6, 129) F 
JX = 30 
JY = JX - 1 
CALL CALMU! UY ) 

JX = JY + 1 

IF { JY .EO. 29 ) WRITE (6, 125) 

C ALL ABSORB ( AN , B , WA , AL , AK , UHAT , GAMMA) 

WRITE (6, 130) B , WA , UHAT . AL 
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DAMP = A1MAG !W(JX) ) 

DSTAR = DAMP’ SOS/RAD 1 US 
WRITE (6, 107) DAMP , DSTAR 

CALL DEFFEC ( W(UX), 8, WA.AL, GAMMA, AK, RADIUS, SOS.RLMDA) 

WRITE (6, 123) RAD , RADIUS 

IF ( MATRIX .NE. 3HYES ) GO TO 1 

DO 6 J = 1 , JY 

KS = 1 

KF=3 

WRITE (6, 106) U, KS, KF 
D06 L=1 ,30 

6 WRITE (6, 104) L, ( MUU.L.J), I = KS, KF ) 

1 CONTINUE 
GO TO 8 
9 CALL EXIT 

100 FORMAT!' 1 THIS IS ENGINE NUMBER \A10) 

101 FORMAT CO THE PRIMARY MODE ASSUMED IS LHAT = ',12,' MHAT = ',12,' 
INHAT = ', 12) 

102 FORMAT CO BETAN = \2G21.14) 

103 FORMAT C+',70X,' N = '.2G21 . 14, 

104 FORMAT C *, 13,' * , 1 0G13.6) 

105 FORMAT CO BETAWl = '.2G21.14) 

106 FORMAT ( *0 N J EQUAL ',12,' L EQUAL ',11,' TO ’.12) 

107 FORMAT ( '0 THE DAMPING RATE IS 'G21.14.28X' THE LOG DAMPING RATE IS 
1 *G2T .14* INVERSE SECONDS') 

108 FORMAT t' + ',70X,' K = \2G2U4) 

109 FORMAT CO OMEGA! 1} = *2021.14) 

•112 FORMAT ('0 BET AW ID = * ,2G21 .14) 

113 FORMAT ('0 W/WO = '2G21.14) 

114 FORMAT !'+',70X,' A = ',2021.14) 

115 FORMAT ( *+* ,70X, ' Th£ LINER STARTS AT 'F6.4' AND ENDS AT *F6.4) 

116 FORMAT CO THE MACH NUMBER IS \G2U4) 

117 FORMAT C + \70X,' THE RATIO OF SPECIFIC HEATS IS * ,G21 . 14) 

118 FORMAT CO THE FREQUENCY IS ', 2G21.14) 

119 FORMAT ('+• ,70X, ' THE LENGTH OF THE COMBUSTOR IS \G2I.14) 

120 FORMAT (*+' ,70X, ' G = *, 2G21.14) 

1 21 FORMAT CO BETA C = ’, 2G2U4) 

122 FORMAT (’0 IMPEDANCE OF THE INJECTOR IS '2F18.14) 

123 FORMAT ( *0 THE RADIUS IS -'F10.5* INCHES OR 'F10.5' FEET') 

1 25 FORMAT ( ' 1 ' , 1 35 !1 HI ) , // , * THE PROBLEM D I D NOT CONVERGE * , //, 1 

1 35 ! 1 H* ) ) 

126 FORMAT CO OMEGA WIGGLEY » '2G21 . 14, 1 1X, ' W/WO WIGGLEY = *2021.14) 
129 FORMAT C+'70X* W/WO = *2021.14) 
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130 FORMAT ! *0 THE BACKING DISTANCE IS * G21.14.24X,* THE'WIDTH OF Tl-£ 
1APERATURE IS 'G21.14//* THE PEAK VELOCITY IS *2G21 . 14.6X, * T*£ LE 
1NGTH OF THE APERATURE IS * G21.14) 

131 FORMAT ! *0 THE DIMENSIONAL FREQUENCY IS 'G21.14J 
200 FORMAT ( A 1 0 , 6F 10.0 /2G2 1 . 14, 212, A4.3F5. 0,215) 

END 

SUBROUTINE ABSORB! A, B,WA, LA, Z.UHAT, GAMMA) 

REAL LA , leff 
COMPLEX A, F , Z , UHAT 

COMMON / BLKI / F.RHOCHAM.RHOVOL.PRESURE.RADIUS.SOS 
W = 1 .841 1 8378 * REAL ( F ) 

LEFF = WA / 2. /W/W/B/B 

LA = LEFF - 0.85 ' WA * ( 1 .-0.7'SQRT (WA/2./B) ) 

IF ( LA .GE. 0.001 / RADIUS 12,1 

1 LA = 0.001 / RADIUS 
XA = w * LEFF 

XC = WA / 2. / W / B / B 
Z1 = XA - XC 

ZR = SQRT ( -Z I’ Z1 +SQRT (ZTZI’ZrZ I +0. 04/GAMMA/GAMMA) l/SQRT (2.0) 

Z = CMPLX ( ZR , Zi ) 

GO TO 3 

2 Z = CMPLX ( SQRT (0.1 /GAMMA) , 0.0 ) 

3 UHAT = CMPLX ( 0.1 / CABS(Z) , 0.0 ) 

F = CMPLX ( W , 0.0 ) 

RETURN 

END 

SUBROUTINE DEFFEC ( W,B, WA, LA, GAMMA, ZI .RADIUS, SOS, RLMDA) 
real la.leff.kr.ki 
COMPLEX W , ZI , Z2, Z3 

KR ( ZI , GAMMA ) = SQRT (-Z! *Z I +SQRT (ZI *Z I ‘Zl’ZI+O. 04/GAMMA/GAMMA) 
1) /SQRT (2.0) 

K I ( W1, LEFF , WA , B ) = W 1 * LEFF - WA/2 . /W 1 /B/B 
LEFF = LA + 0.85’WA* (1.0-0. 7*SQRT(WA/B5) 

WR = REAL (W) 

Wl = A 1MAG (W) 

Z! * K1 (WR* 1.5, LEFF, WA,B) 

Z2 = CMPLX ( KR(Z1, GAMMA) ,Z1) 

ZI = K 1 (WR/1.5,LEFF,WA.B) 

Z3 = CMPLX ( KR (ZI , GAMMA) , Z!) 

DAMP = Wl*3.0/d.0+(CABS(Z2)+CABS(Z3))/CABS(ZD) 

DSTAR = DAMP* SOS/RAD I US 
WRITE (6, 100) DAMP , DSTAR 

100 FORMAT CO Th€ EFFECTIVE DAMPING RATE IS ’G21.14.18X* THE EFFECTIVE 
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1 L, D. R. IS ‘G21.14* INVERSE SECONDS') 

RETURN 

END 

SUBROUTINE CUBIC IP , Q , R , B. DELSTAR) 
AX=(3.’Q-P*P)/3. 

BX-- (2. ’P’P‘P-9. *P*Q+27. *R)/27. 

AX 1 = AX V AX * AX / 27. 

8X1= BX * BX / 4. 

D1S = AX1 + BX1 

IF ( ABS ( DIS / AX1 ) .LE. 1.E-07 ) DIS = 0.0 
IF ( DIS .LT. 0.0 ) GO TO 1 
IF (DIS) 1 ,2,2 

1 PH I = ( AC OS (-BX/2 . /SORT ( - AX * AX * AX/27 . ) ) ) /3 . 
XC0EFF=2. ’SORT ( -AX/3.) 

CON=0. 017453292519943 
X1=XC0EFF*C0S(PH|)-P/3. 

X2= XC OEFF * C OS ( PH 1 + 1 2 0 . ' C ON ) -P/3 . 

X3=XC0EFF*C0S (PH 1+240. *C0N) -P/3. 

IF (XI .LT.O.O. OR .XI . G T . B )X1=-1000. 

1F(X2. LT.O.O. OR. X2.GT. B )X2=-1000. 

IF(X3. LT.O.O. OR. X3.GT. B )X3»-1000. 
DELSTAR=AMAX1 (XI ,X2,X3) 

GO TO 4 
2EXP=1./3. 

BX2=-BX/2. 

C1=BX2+SQRT (DIS) 

C2=BX2-SQRT(DIS) 

ICODE=0 

JCODE=0 

IF (Cl .LT.O.O) 1C0DE=1 
IF (C2.LT.O.O)UCODE=1 
AXE=(ABS(CD)’*EXP 
BXE= (ABS (C2) ) * 

IF ( ICODE) 7,8,7 

7 AXE=-AXE 

8 IF (JC0DEJ9, 10,9 

9 BXE=-BXE 

10 IF (DIS) 1.1,11,3 

11 X1=A/E+BXE-P/3. 

•/ 2 =-x 1 /2 . -P/2 . 

if(xi.lt!o.o.or.xi.gt. B )X1=-1000. 

IF(X2. LT.O.O. OR. X2.GT. B )X2=-1000. 
DELSTAR=AMAX1 (XI ,X2) 
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GO TO 4 

3 DELST AR=AXE+BXE~P/3 . 

4 IF (DELST AR.LE. 0.0. OR. DELSTAR.GE. B ) GO TO 5 
RETURN 

5 WRITE (6,6) 

CALL EXIT 

6 FORMAT { 1H1 , * THIS VALUE IS OUT OF RANGE.*) 

END 

FUNCTION BESPR I M (M, L) 

DIMENSION A { 1 0 ,5) 

C*”* THESE ARE THE ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF OR 
C”“ SET EQUAL TO ZERO 

A(1 ,1) = 0.00000000 
A (2 , 1) = 3.83170597 
A (3 ,1) = 7.01558667 
A (4 ,1) = 10.17346814 
A (5 ,1) = 13.32369194 
A (6 ,1) = 16.47063005 
A (7 ,1) = 19.61585851 
A (8 ,1) = 22.76008438 
A (9 ,1) = 25.90367209 
A ( 1 0 , 1 ) = 29.04682853 

C*’** THESE ARE THE ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF OR 
C*'** SET EQUAL TO ZERO 

A ( 1 ,2) = 1.84118378 
A (2 ,2) = 5.33144277 
A (3 ,2) = 8.53631637 
A (4 ,2) = 11.70600490 
A (5 ,2) = 14.86358863 
A (6 ,2) = 18.01552786 
A (7 ,2) = 21.16436986 
A (8 ,2) = 24.31132686 
A (9 ,2) = 27.45705057 
A (10,2) = 30.60192297 

“** THESE ARE THE ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF OR 
**•• SET EQUAL TO ZERO 
Ad ,3) = 3.05423693 
A(2 .3) = 6. "’0613319 
A (3 ,3) = 9.96946782 
A (4 ,3) = 13.17037086 
A (5 ,3) = 16.34752232 
A (6 ,3) = 19.51291278 
A (7 ,3) = 22.67158177 
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AiB ,3) = 25.82bQ3714 
A (9 ,3) = 28.97767277 
A (10,3) = 32.12732702 

C”" THESE ARE THE ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF OR 

C*“* SET EQUAL TO ZERO 

All ,4) = 4.20118894 
A (2 ,4) = 8.01523660 
A (3 ,4) * 11.34592431 
A (4 ,4) = 14.58584829 
A (5 ,4) = 17.78874787 
A (6 ,4) * 20.97247694 
A (7 ,4) = 24.14489743 
A (8 ,4) « 27.31005793 
A (9 ,4) 8 30.47026881 
A(10, 4) = 33.62694918 

C“" THESE are the ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF OR 

C'“* SET EQUAL TO ZERO 

A(1 ,5) = 5.31755313 
A (2 ,5) = 9.28239629 
A (3 ,5) = 12.68190844 
A (4 ,5) = 15.96410704 
A (5 ,5) 8 19.19602880 
A (6 ,5) = 22.40103227 
A (7 ,5) = 25.58975968 
A (8 ,5) = 28.76783622 
A (9 ,5) = 31.93853934 
A ( 10.5) = 35.10391668 
1 BESPRIM = A (L, M) 

RETURN 
ENTRY BESSEL 

C*'“ THESE are the BESSEL NUMBERS OF ORDER ZERO FOR THE ZEROS OF THE B 

C"** FUNCTION 

A(1 ,1) = 1.00000000 
A 12 ,D = -0.4027588095 
A (3 ,1) = 0.301128303 
A (4 ,1) = -0.249704877 
A (5 ,1) = 0.218359407 
A (6 ,1) 8 -0.19645371 
A (7 ,1) = 0.180063375 
A (8 ,1) = -0.167184600 
A (9 ,1) = 0.156724985 
A ( 1 0 , 1 ) = -0.148011108 

c * " * these are the bessel numbers or order one for the zeros of the b 
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C*“* FUNCTION 

A(1 ,2) = 0.5818649368 
A (2 ,2) = -0.346)258542 
A (3 ,2) = 0.273298113! 

A (4 ,2) = -0.233304416 
A (5 .2) = 0.207012651 
A (6 ,2) = -0.188017488 
A (7 ,2) = 0.173459050 
A (8 ,2) = -0.161838211 
A (9 ,2) = 0.152282069 
A (10,2) = -0.144242905 

C"" THESE are THE BESSEL NUMBERS OF ORDER TWO FOR THE ZEROS OF THE B 
C**“ FUNCTION 

A(1 ,3) = 0.4864961885 
A (2 ,3) = -0.3135233099 
A (3 ,3) = 0.2547441235 
A (4 ,3) = -0.220881581 
A (5 ,3) = 0.197937434 
A (6 ,3) = -0.181010000 
A (7 ,3) = 0.167835534 
A (8 ,3) = -0.157195167 
A (9 ,3) = 0.148363778 
A (10, 3) = -0.140878333 

c"** these ARE the BESSEL NUMBERS OF ORDER THREE FOR THE ZEROS OF THE B 
C'**’ FUNCTION 

Ad ,4) = 0.4343942763 
A (2 ,4) = -0.2911584413 
A (3 ,4) = 0.240738175 
A (4 ,4) = -0.210965204 
A (5 ,4) = 0.190419022 
A (6 ,4) = -0.175048405 
A (7 ,4) = 0.162954965 
A (8 ,4) = -0.153102409 
A (9 ,4) = 0.144866574 
A ( 10,4) = -0.137844513 

C**” THESE are THE BESSEL NJMBERS OF ORDER FOUR FOR THE ZEROS OF THE B 
C"“ FUNCTION 

Ad ,5) = 0.3996514545 
A (2 ,5) = -0.2743809949 
A (3 ,5) = 0.229590468 
A (4 ,5) = -0.202763849 
A (5 ,5) = 0.184029896 
A (6 ,5) = -0.169878516 
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A (7 ,5) = 0.158655372 
A (8 ,5) = -0.149451156 
A (9 ,5) = 0.141714307 
A ( 1 0,5) = -0.135086328 
GO TO 1 
ENO 

SUBROUTINE CALCA4( A4, A5 ) 

OmEX BETA1 , BETA2, BETAN, BETAC, BETA!, BETAWI, OMEGA, C, AK1, 

1 A K2, PS I 2, SOME, SOME 1 
COMPLEX A4, A5, PS1, DUMMY, G1, G2, 1 

COMMON / BLKA / BETA1 , BETA2, AMACH, GAMMA, OMEGA, ALENGTH, C, AK1 
1, AK2, BETAN, RLK)A, BETA I, BETAC, BETAWI, 1, LHAT, MHAT, NHAT, 

1 PSI2 ,PI,X2,X1,S0ME,S0ME1 

A4 = AIC1/EPSIZN(LHAT, MHAT, NHAT, ALENGTH, RLMDA1/PSI (DUMMY) 

A5=BETAC * (G 1 (X2.MHAT) -G 1 ( X 1 , NHAT 5) /EPS I KAT (LHAT , MHAT , NHAT , ALENGTH, 
1RLMDA1/PS1 (DUMMY) 

RETURN 

END 

FLECTION EPSIKAT (L,M,N, ALENGTH, RLMDA) 

DUM1 = ALENGTH 

IF ( N .NE. 0) DUM1 = DUM1/2. 

IF C M .EQ. 0) GOTO 1 

EPSIKAT ■ (0.5-M *M / (2. *RLK)A*RLMDA) ) *DUM1 
RETURN 

1 EPS KAT = DUMI/2. 

RETURN 

ENTRY EPS1ZN 

IE( N .EQ. 0 ) GO TO 2 

EPSIKAT* ALENGTH / 2. 

RETURN 

2 EPSIKAT* ALENGTH 
RETURN 

ENTRY ETA 

EPSIKAT * (N-I)'(N-I) ' 3.1415926 *3.1415926 / (ALENGTH * ALENGTH) 
1 +RLMDA ’ RLMDA 
RETURN 
DC 

F'Jf JC T I ON EPS 1 R2N ( L , M , N , ALENGTH, RLMDA) 

COMPLEX FN0RM,DW2 
DUM1 « 1. 

IE(M.EQ.O) D'JMI - .5 

DUM2-D JORM ( L , M f N , RLK) A , ALDJG TH ) ' ' 2 

EPS I RZN-DUM 1 ' REAL (DUM2 ) /3 . 1 4 1 5926 
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RETURN 

END 

COMPLEX FUNCTION FNORM ( L , M , N , RLMD A , ALENG TH ) 

DUM1 = 3.1415926 

IFi M .EQ. 0 ) DUM1 = 6.2831852 

DUM 1 =DUM 1 * EPS I KAT (L , M , N , ALENG TH , RLMD A 5 * BESSEL (M+1 ,L) * *2 
FNORM=CMPLX (DUM1 ,0.0) 

FNORM= CSQRT (FNORM) 

RETURN 

EM) 

COMPLEX FUNCTION G2(N) 

COMPLEX BETA1 , BETA2, BETAN, BETAC, BETA!, BETAWI , OMEGA, C, AK1 , 

1 AK2, PS12, MU, W, F, I, SOME, SOME1 
COMPLEX DUM1, GDUM, BETA, GBUM, DUM 

COMMON / BLKA / BETA 1 , BETA2, AMACH, GAMMA, OMEGA, ALENGTH, C, AK1 
1, AK2, BETAN, RLMDA, BETA!, BETAC, BETAWI, I, LHAT, MHAT, NHAT, 

1 PS 12 ,P1,X2,X1,S0ME,S0ME1 

GBUM (BETA, DUM1.N, AMACH, OMEGA, ALENGTH) = - (AMACH ’BETA 'BETA+2. ‘OMEGA 
1 ‘BETA) / ( ( N ‘ 3.1415926 / ALENGTH ) “2 - BETA* *2 ) *(CEXP(DUMI* 

1 ALENGTH) * DUM1 * M.)**(N+2) - DUM1 ) 

DUM (DUM1, ALENGTH, NHAT) = CEXP (DUM1 ‘ALENGTH) * (-1 .)** (NHAT+2) -1 . 

BETA = BETA1 

DUM1 = CMPLX(-AIMAG(BETA1),REAL(BETAD) 

G2 = GBUM ( BET A , DUM 1 , N , AMACH , O^GA , ALENGTH) 

BETA = BETA2 

DUM1 = CMPLX(-AIMAG(BETA2),REAL(BETA2)) 

G2 = G2 + C * GBUM (BETA, DUM1 ,N, AMACH, OMEGA, ALENGTH) 

IF ( N .EQ. 0 .AM). CABS ( BETA2) .LE. l.E-10 ) RETURN 

RETURN 

ENTRY PS I 

DUM1 = CMPLX( -AIMAG (BETA1 ) , REAL (BETA 1 ) ) 

G2 = DUM1*DUM(DUM1, ALENGTH, NHAT) / (WAT* *2*3. 1415926* ’2/ALENGTH* *2 
1 -BETA 1* *2) 

IF ( NHAT .EQ. 0 .AM). CABS ( BETA2) .LE. l.E-10 ) GO TO 3 
DUM1 = CMPLX ( -AIMAG ( BET A2), REAL (BETA2) ) 

G 2 = ( G2 +C’DUM1’DUM(DUM1, ALENGTH, NHAT)/(WAT**2*3. 1415926* *2/ 

1 ^ ALENGTH* * 2-BETA2 ’ * 2 ) ) /EPS 1 ZN (LHAT , MHAT , NHAT , ALENGTH , RLMDA) 

RETURN 

3 G2 = ( G2 + C ’ ALENGTH ) / EPS I ZN ( LHAT , MHAT , NHAT , ALENGTH, RLMDA ) 
RETURN 
DO 

SUBROUTINE CALMU(JP) 

COMPLEX aa,ab,ac,ad,ae,af,ag,ah,ai,aj,ak,al,am,an,ao,aiii 
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COMPLEX A1 , A2,A3,A4, A5,A6,A7,A8, A9.SUM1 ,SUM2,SUM3,SUM4,G1 ,G2,PS1 
COMPLEX BETA1, BETA2, BETAN, BETAC, BETA I, BETAWl , OMEGA, C, ,AK1, 

1 AK2, PS12, MU, W, E, |, ENORM, D1,D2,AK1J,AK2J,SSG1,A11,S0MEiS0ME1 
COMMON / BL<A / BETA1 , BETA2, AMACH, GAMMA, OMEGA, ALENGTH, C, AK1 
1, AK2, BETAN, RLMDA, BETA!, BETAC, BETAWl, 1, LHAT, MHAT, NHAT, 

1 PS12 , P I , X2 , X 1 , SOME , SOME 1 
COMMON / BLKB / MU( 3,30,29) , W(30) 

COMMON Cl 1 (30,30) , 00(30,30) , ERROR 1 , ERR0R2 , LDEX , M)EX 
DKEXP(I’ALENGTH’BETAI) 

D2=CEXP ( I *ALENGTH*BETA2) *C 
A 1 =W ( 1 ) **2 
A4=0MEGA* *2 

A2 = BETAC ’BESSEL (MHAT+1, LHAT) /PS 1 2 /ENORM (LHAT , MHAT , NHAT , RLIDA , A 
1 LENGTH) 

• AK2 J= 1 * GAMMA ’ (Wl 1 ) * (D1+D2)+AMACH* (BETA1 *D1+BETA2'D2) ) 

AK1J=] ’GAMMA’ (W( 1 ) * ( 1 . +C) +AMACH* (BETA1 +C*BETA2) ) 

DO 1 L = 1 ,LDEX 

RLMDB = BESPR I M (MHAT+1, L) 

A3 = BESSEL (MHAT+1 ,L) ’CMPLXd .0,0.0) 

DO 1 N = 1 ,NDEX 
N1 = N-1 

A 1 = ( A4 - ETA (L, MHAT, N, ALENGTH, RLMDB)) 

AH= ( A 1 - ETA (L, MHAT, N, ALENGTH, RLMDB)) 

IF (L.EQ.LHAT. AMD.N1 .EQ.NHAT) GOTO 1 

AO=OMEGA 

0MEGA=W ( 1 ) 

AF=G1 (X2.N1) - G1 (XI ,N1 ) 

AG=G2(N1) 

0MEGA=A0 

AN=CMPLX(0. 0,0.0) 

IE (L.NE.LHAT) GOTO 31 

AN=BETAI’A<1J + BETAN*AK2J’(-1.-)”N1 + AMACH’AG 
AN=AN* (A4 -AD /AH 

AN=AN + 1*(W(1) - OMEGA)* (GAMMA’ (BETAI* (1. + C) + BETAN* ID 1 +D2) 

1 ’ (-1 . ) ”N1 ) + AMACH’SSGI (ALENGTH, ND) 

AN= AN/PS 1 2 /ENORM (LHAT , MHAT , NHAT , RLMDA , ALENGTH) /EPS I ZN (LHAT , MHAT , 
INI, ALENGTH, RLMDA) /A I 
31 CONTIS 

MU(L,N,1) = (A2/A3)’AE/A1/EPSKAT(L, MHAT, N1, ALENGTH, RLMDB) + AN 
1 CONT 1MJE 
NHA1 = NHAT+1 

KJ (LHAT, hJHA 1,1) =CMPLX ( 0.0 ,0.0 ) 

DO 2 NX = 1 ,NDEX 
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N = NX-1 

DO 2 NY = 1 ,NDEX 
NP = NY - 1 
IF (N.EQ.NP)3,4 

3 1F(NP.EQ. 0)5,6 

4 Cl 0 (NX, NY) =GAMMA/2./ IN* '2-NP* '2) ' ALENGTH /PI* UN+1P) * (S1N( 

1 (N— NP) *PI ’X2/ALENCTH) -SIN! (N-NP) ‘PI 'XI /ALENGTH) ) + (N-NP) * (SIN( (N+NP 
1 ) *P 1 * X2/ALENGTH) -S 1 N ( (N+NP) ’ P I* X 1 /ALENGTH) 5 ) 

C 11 (NX, NY) = GAMMA / 2. / ( N”2 - NP* *2) * AMACH’NP* ( (NP+N) * (COS 
1 ( (NP-N) ’PI * XI /ALENGTH) -COS ( (NP-N) ’PI ’X2/ALENGTH) ) + (NP-N) ’ (COS ( (NP+ 
IN) *P I ’XI /ALENGTH) -COS ( (NP+N) *P1’X2/ALENGTH))) 

GO TO 2 

5 C 1 0 (NX, NY) =GAMMA* (X2-X1) 

Cl 1 (NX, NY) =0.0 

GO TO 2 

6 CIO (NX, NY) =GAMMA* (ALENGTH/ (4. *N’PI ) ’ (SIN (2. ’N’PI *X2/ ALENGTH) 
1-S1N(2.’N’PI’X1/ALENGTH)) + (X2-XD/2.0) 

Cl 1 (NX, NY) = - GAMMA‘AMACH/2. ‘(SIN(N*P|’X2/ALENGTH)”2-S1N(N*PI*X1 
1 /ALENGTH) ”2) 

2 CONTINUE 
CALL W1 (2) 

IF ( JP ,EQ. 1 ) RETURN 
DO 7 J = 2, UP 
A1 = W(U) - W(U-I) 

IF ( ABS ( REAL ( A 1 ) /RE AL ( W ( J- 1 ) ) ) . LE . ERROR 1 . AM) . ABS ( A I MAG ( A 1 ) / A I MAG ( W 
1 (J-1 ) ) ) .LE.ERR0R2) GO TO 13 
A1=W( J )”2 
AU=0MEGA”2 

A2=BETAC ’BESSEL (MHAT+1 ,LHAT)/PSI2 /FNORM (LHAT , MHAT , NHAT , RLMDA , ALE 
1NGTH) 

A3 = BETAI ’GAMMA* l’W( U ) 

A4 = BETAN’GAMMA* I *W(J ) 

A5 = 2. * I *W(J ) ’AMACH 

AK1J=1 ’GAMMA’ (W(J ) * ( 1 ,+C) +AMACH* (BETA 1 +C ’BET A2) ) 

AK2U=1 ’GAMMA* (W(J ) * (D1+D2)+AMACH* (BETA1 *D 1 +BETA2 *D2 ) ) 

A6 = (AMACH’PI) ”2/2. /ALENGTH’CMPLXd. 0,0.0) 

DO 8 L = 1 , LDEX 

RLMDB = BESPR I M (MHAT + 1 , L) 

A7 = CMPLX (BESSEL (MHAT+1, L) , 0.0 ) 

DO 8 N = 1 ,NDEX 
fJ1 = N-1 

IF ( L.EQ.LHAT .AM). N1 .EQ. WAT) GO TO 8 
S'JMI = S'JM2 = SUM3 = SUM4 = CMPLX (0.0, 0. 0) 
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DO 9 LP = 1.LDEX 

A8 = BESSEL (MHAT+1 ,LP) 'CMPLXU .0,0.0) 

DO 10 NP = 1.NDEX 
NP I =NP- 1 

10 SUM3=SUM3+MU(LP,NP,J-D *(C10(N,NP) TWIJ )+C1 1 (N,NP) ) *A8 
9 CONT 1 NUE 

DO 11 NP = 1.NDEX 
A8=MU( L , NP , U- 1 ) 

SUM1 = SUM1+A8 
SUM2=SUM2+A8* (-1 . ) * ' (NP+1 ) 

IF (N.EQ.NP) 11,12 

1 2 SUM4=SUM4+A8* (M>- 1 ) “2* U .- {- 1 . ) * * (NP+M) ) / ( (NP- 1 ) ' *2-M1 * ‘2) 

11 CONTINUE 

A8 = EPSIKAT (L ,MHAT,N-1,ALENGTH,RLMDB)*CMPLX(1.0,0.0) 

A9=EPS 1 ZN (LHAT , MHAT , N- 1 , ALENGTH, RLMDA) ‘CMPLX ( 1 .0, 0. 0) 

A IOMEGA 
OMEGA=W(J ) 

AJ=G2(N1) 

AA = G1 (X2.ND-G1 (XI, N1) 

OMEGA=A I 

AB = A2'AA/A8/A7 

AC = A3* SUM 1 /A9 

AD = A4*M.0)*'(N+1)'SUM2/A9 

AE = BETAC/A8/A7‘SUM3 

AF = -A5/A9*SUM4 

AG = -A6/A9*N1 * *2*MU(L,N f U-1 ) 

AI= A1 1 -ETA (L,M,N, ALENGTH, BESPRJM (MHAT+1 ,L) ) 

AH=A1 -ETA (L.M.N, ALENGTH, BESPRm(MHAT+1,L)i 
AM= (AB+AC+AD+AE+AF+AG) /AH 
AN= CMPLX(0. 0,0.0) 

IF(L.NE.LHAT) GOTO 400 

AN=(BETA!*AK1J + BETAN*A<2J* (-1 . ) * ’ (N+1 )+ AMACH* A J) /PS 1 2/FNORM ( 

1 LHAT , MHAT , NHAT , RLMDA , ALENGTH) /A9/ ( AH* A I ) 

AN=AN* (All -A1) 

AO=SSG 1 ( ALENGTH, N1) 

AN=AN+(W(J ) -OMEGA) * ( I ‘GAMMA ’BETA I * (1 .+C) + I ‘GAMMA’BETAN* (D1+D2) 
1* ( - 1 . ) * *N1 + 2 . * I ’ AMACH* AO) /PS 1 2/FNORM (LHAT , MHAT , NHAT , RLMDA , ALENGTH 
2)/A9/Al 
400 CONTINUE 
8 MU(L,N, J) = AM+AN 
fjl = NHAT+1 

MI!LHAT,fJ1,J)=CMPLX(0.0,0.0) 

7 CALL W1 (J+1) 
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RETURN 

13 JP = J - 1 
RETURN 
END 

SlBROUTlte W1 (J) 

COMPLEX SUM1 ,SUM2,SUM3,SUM4,A1 ,WlETA,G1 ,G2,PS1 ,FNQRM, A.DUM1 ,DUM2 
COMPLEX BET A 1 , BETA2, BETAN, BETAC, BETA!, BETAWi, OMEGA, C, AK1, 

1 AK2, PS12, MU, W, F, 1 , BET A, SOME, SOME 1.SUM5, ODUM, E,G 
COMMON / BLKA / BETA 1 , BETA2, AMACH, GAMMA, OMEGA, ALENGTH, C, AK1 
1, AK2, BETAN, RLMDA, BETA!, BETAC, BETAWI, I, LHAT, MHAT, NHAT , 

1 PSI2 , P I , X2 , X 1 , SOME , SOME 1 
COMMON / BLKB / MU( 3,30,29) , W1305 

COMMON Cl 1 (30,30) , CIO (30,30) , ERROR! , ERR0R2 , LDEx , NDEX 
SUM! = SUM2 = SUM3 = SUM4 = SUM5=CMPLX (0. 0, 0. 0) 

DO 1 L = 1 , LDEX 

AKMPLX (BESSEL (MHAT+1.L) , 0.0) 

DO 1 N = 1 ,NDEX 

E=MU(L,N, J-1 ) ’A! 'CIO (NHAT+1 ,N) 

SUM5=E + SUM5 

SUM4=E'|*W(J-1) + MU(L,N, J-1 ) *A1 *C1 1 (NHAT + 1,N) + SUM4 

1 CONTINUE 

DO 2 N1 = 1.NDEX 
N = Ni - 1 

IF (N.EQ.NHAT) GO TO 2 
A! = MU (LHAT, NI, J-1) 

SUM1 = SUM1+A1 

SUM2 = SUM2 + A1'M.)"N 

SUM3 = SUM3+A1 *N* *2/ (N"2“NHAT* '2) * (1 . - (-1 . ) " (N+NHAT) ) 

2 CONTINUE 
BETA=OMEGA 
OMEGA=W(J-1) 

W I ETA=FNORM (LHAT , MHAT ,INHAT, RLMDA , ALENGTH) * ( I 'OMEGA/EPS I ZN (LHAT , MHA 
IT, WAT, ALENGTH, RLK5A) ' (GAMMA* (BETA! 'SUM1+BETAN* (-1 . ) " (NHAT+2) 'SUM 
12) -2. 'AMACH*SUM3)+BESSEL(MHAT+1 ,LHAT) 'BETAC/EPS I RZN (LHAT, MHAT, NHAT 
1 .ALENGTH, RLK)A)*SUM4) 

OMEGA= BETA - 

W1ETA = WIETA+ (W( 1 ) "2 + ( W ( J— 1 ) — W ( 1 ) J'SOME) 

DDUM=WIETA 
F <SQRT (WiETA) 

E=F 

IF (REAL (F ) .LT.0.0) F=-F 

G= (BESSEL (MHAT + 1 , LHAT ) 'BET AC/EPS ! RZN (LHAT , MHAT , NHAT , ALENGTH, RLK)A 
1 ) 'S'jra+ { GAMMA ' BET A 1 ' SUM 1+G AMMA ' BETAN' (-1 . ) ' *NHAT-AMACH'2. *SUM3) /E P 



C-32 


1 S l ZN (LHAT , HHAT , NHAT , ALENGTH, RLMD A ) ) *FNQRM (LHAT, MHAT, NHAT, RLM)A, ALE 
1NGTH) 

G = G * i + SOME 
NN = 0 
20 CONTI NX 

W1ETA = DDUM+ (F-WIJ-1 ) ) *G 
F = CSQRT (U1ETA) 

IF { REAL (F) .IT. 0.0 ) F = -F 
IF ( CABS (F-E) ,LT. 1 .E-07 .OR.NN.GT. 105 12, 13 
13 E = F 
NN = NN+1 
GO TO 20 
12 W(J) = F 

WRITE(6, 100) J, W(J) 

F = F / 1.84118378 
WRITE (6, 101 ) F 
RETURN 

100 FORMAT (*0 OMEGA (M2,*) = \2G21.14) 

101 FORMAT (*+‘70X* W/W0 = ’2G21.14) 

EM) 

COMPLEX FUNCTION G1 ( X , N ) 

COMPLEX BETA1. 8ETA2, BETAN, BETAC, BETA I , BETAWI , OMEGA, C, AK1 , 
1A<2,PS12,MU,W,F,I,BETA,S0ME ( S0ME1 
COMPLEX DUM1 , GOUM, GBUM, DUM 

COMMON / BL<A / BETA1 , BETA2, AMACH, GAMMA, OMEGA, ALENGTH, C, AK1 
1, AK2, BETAN, RLMDA, BETA), BETAC, BETAWI, I, LHAT, MHAT, NHAT, 

1 PSI2 ,PI,X2,X1,SOME,SOME1 
GOUMIA.OUMI .BETA ,X)=1 ./ (A* ’2-BETA* '2 
1 ) * (CEXP (DUM1 * X) * (DUM1 *COS(A*X) +A*S]N(A'X) ) -DUM1 ) 

IF ( X .EQ. 0.0 ) GO TO 1 
A = N*3. 14 15926/ ALENGTH 

G1=GDUM(A, l’BETA1,BETAl,X)' (OMEGA+AMACH'BETAI) 

IF ( N .EQ. 0 .AM). CABS l BETA2) .LE. 1.E-10 5 GO TO 2 
G1 = GAMMA * (G1 + C * GDUMIA, I'BETAI ,BETA1 ,X) * (0MEGA+AMACH*BETA2 
1)5*1 
RETURN 

1 G1 = CMPLX (0.0, 0.0) 

RETURN 

2 BETA = C * OMEGA * X 

01 = GAMMA * ( G1 + BETA) ’] 

RET'JRN 
ENTRY SGI 

IF l X .EQ. 0.0 ) GO TO 11 
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A = N*3. 1415926/ALENGTH 
G1=G0UM(A, I’BETAl f BETA 1 , X) 

IF ( N .EQ. 0 .AND. CABS { BETA2) .LE. 1.E-10 ) GO TO 21 
G 1 = G1 + C * GDUM(A, 1 ‘BETA 1, BETA), X) 

RETURN 

11 G1=CMPLX(0. 0,0.0) 

RETURN 
21 BETA=C'X 
G 1 =G 1 +BETA 
RETURN 
ENTRY SSG1 

IF ( X .EQ. 0.0 ) GO TO 111 

A = N* 3. 1415926/ALENGTH 

G1=GDUM(A, I *BETA1 , BET A 1 , X ) ' I 'BETA1 

1R N .EQ. 0 .AND. CABS! BETA2) .LE. l.E-10 ) RETURN 

G1 = G1 + C ' GDUM(A,i*BETAl,BETAl,X) , l'BETA2 

RETURN 

111 GKMPLXtO. 0,0.0) 

RETURN 

END 


THE FOLLOWING IS AN EXAMPLE OF THE PUNCHED CARD 1M>UT. 

ONE 1.984159 0.265865 9.07012 5181.661961.146 0.982 

+7.613141 23437926E-02+6 . 37493806 1 92886E-02 1 0 NOO. 1000.0100.010 3 3 C 
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THIS IS ENGINE NUMBER ONE 

THE MACH NUMBER IS .26586500000000 

THE RATIO OF SPECIFIC HEATS IS 1.1460000000000 

THE FREQUENCY IS 1.8080424719600 3.23641657590519E-14 

THE LENGTH OF THE COMBUSTOR IS 1.9841590000000 

THE PRIMARY MODE ASSUMED IS LHAT = 1 MHAT = 1 NHAT = 0 

THE LINER STARTS AT .0965 AND ENDS AT .1035 

A = 1.0000000000000 0. 

BET AN = 1.16979672024435E-02 0. 

G = .91660000000000 0. 

BET AC = 2.9539809563369 0. 

t' = ?Q c i , SQRflQ f i6 5 ? ; ?7fl 0 

BETAWl’* 7.61314123437886E--02 6. 374938061 93 157E-02 

N = .58624670209981 -.23978101901084 1 

BETAW1D = 7.61314123437034E-02 6. 374938061 920 15E-02 

N = .58624670210013 -.23978101901041 I 

IMPEDANCE OF THE INJECTOR IS .07613141234379 .06374938061929 

N = .58624670209980 -.23978101901073 I 

OMEGA W1GGLEY = 1 . 80804247 19600 3 . 2364 1 6575905 1 9E- 1 4 

W/WO W I GGLEY = .981 99999999997 1 . 75779 1161 89323E- 1 4 

THE BACKING DISTANCE IS .10000000000000 

THE WIDTH OF THE APERATURE IS 6.92336549877837E-03 

THE LENGTH OF THE APERATURE IS .10077525063986 

THE PEAK VELOCITY IS .33852621759621 0. 

THE DIMENSIONAL FREQUENCY IS 1.8080424719600 
OMEGA ( 1) = 1.8077193093524 1 .91 8972408533 18E-02 

W/WO = .98182448106968 1 .0422492471301 IE-02 

OMEGA! 2) = 1.7957796526966 1 . 9901 7000807443E-02 

W/WO = .97533970926934 1 .08091 871 636758E-02 

OMEGA! 3) = 1.7947392480612 1 .9091 16231 59963E-02 

W/WO = .97477463551259 1 .0368960732424 IE-02 

OMEGA! 4) = 1.7937225329955 1 . 89437574447553E-02 

- W/WO = .97422242824421 1 .028890 08965499E- 02 

Th£ BACKING DISTANCE IS .10000000000000 
THE WIDTH OF THE APERATURE IS 6 . 92336549877837E-03 
THE LENGTH OF THE APERATURE IS . 1008131 1485223 
TF€ PEAK VELOCITY IS .33852621759621 0. 

THE DAMPING RATE IS 1 .89437574447553E-02 

THE LOG DAMPING RATE IS 129.86837748249 INVERSE SECOWS 

Th£ EFFECT l VE D AMP t NG RATE IS 1 . 80389999333238E-02 

THE EFFECTIVE L. D. R. IS 123.66583871121 INVERSE SECON)S 

ThC RADIUS IS 9.07012 INCHES OR -.75584 FEET 
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